Log In Sign Up

Matrix Product Operator Restricted Boltzmann Machines

by   Cong Chen, et al.
The University of Hong Kong
Delft University of Technology

A restricted Boltzmann machine (RBM) learns a probability distribution over its input samples and has numerous uses like dimensionality reduction, classification and generative modeling. Conventional RBMs accept vectorized data that dismisses potentially important structural information in the original tensor (multi-way) input. Matrix-variate and tensor-variate RBMs, named MvRBM and TvRBM, have been proposed but are all restrictive by model construction, which leads to a weak model expression power. This work presents the matrix product operator RBM (MPORBM) that utilizes a tensor network generalization of Mv/TvRBM, preserves input formats in both the visible and hidden layers, and results in higher expressive power. A novel training algorithm integrating contrastive divergence and an alternating optimization procedure is also developed. Numerical experiments compare the MPORBM with the traditional RBM and MvRBM for data classification and image completion and denoising tasks. The expressive power of the MPORBM as a function of the MPO-rank is also investigated.


page 6

page 7


Restricted Boltzmann Machine with Multivalued Hidden Variables: a model suppressing over-fitting

Generalization is one of the most important issues in machine learning p...

Data Mapping for Restricted Boltzmann Machine

Restricted Boltzmann machine (RBM) is interpreted as a data mapping betw...

Can a Hebbian-like learning rule be avoiding the curse of dimensionality in sparse distributed data?

It is generally assumed that the brain uses something akin to sparse dis...

Tensor networks for unsupervised machine learning

Modeling the joint distribution of high-dimensional data is a central ta...

A Support Tensor Train Machine

There has been growing interest in extending traditional vector-based ma...

Dimension of Marginals of Kronecker Product Models

A Kronecker product model is the set of visible marginal probability dis...

Matrix Variate RBM and Its Applications

Restricted Boltzmann Machine (RBM) is an importan- t generative model mo...

I Introduction

Arestricted Boltzmann machine (RBM) [1]

is a probabilistic model that employs a layer of hidden variables to achieve highly expressive marginal distributions. RBMs are an unsupervised learning technique and have been extensively explored and applied in various fields, such as pattern recognition 


, computer vision 

[3] and signal processing [4]

. However, conventional RBMs are designed for vector data and cannot directly deal with matrices and higher-dimensional data, which are common in real-life. For example, grayscale images are second-order tensors (i.e. matrices) while color images or grayscale videos are third-order tensors. The traditional approach to apply an RBM on such high-dimensional data is through vectorization of the data which leads to two drawbacks. First, the spatial information in the original data is lost, thus weakening the model’s capability to represent these structural correlations. Second, the fully connected visible and hidden units may cause overfitting since the intrinsic low-rank property of many real-life data is disregarded.

To address these, we propose a matrix product operator (MPO) restricted Boltzmann machine (MPORBM) where both the visible and hidden layers are in tensor forms. The link between the visible and hidden layers is represented by an MPO, which is essentially a tensor network representation of a matrix. To train an MPORBM, we further describe its customized parameter learning algorithms. The MPORBM is also compared with the standard RBM and matrix-variate RBM [5] for tensor inputs in numerical experiments. This article has the following major contributions:

  1. The MPORBM is proposed for the first time and generalizes existing RBM architectures. The standard RBM, matrix-variate RBM (MvRBM) [5] and tensor-variate RBM (TvRBM) [6] models are all special cases of the MPORBM.

  2. Compared with standard RBMs, the number of parameters in an MPORBM grows only linearly with the order of the tensor instead of exponentially. This alleviates overfitting, which is demonstrated through numerical experiments.

  3. Both the visible and hidden layers are represented in tensor forms and therefore useful structural information in the original tensor data is preserved and utilized.

  4. The graphical “language" of tensor network diagrams [7] is introduced to represent how specific quantities (such as the conditional probabilities) are computed. This way, complicated mathematical expressions are easily understood in a visual manner.

  5. Although the data structures are generalized to tensors, the bipartite graph nature of an RBM still applies, together with the fast sampling and inference properties.

Ii Related work

Real-life data are extensively multiway. Researchers have been motivated to develop corresponding multiway RBMs. For example, [8] proposed a factored conditional RBM for modeling motion style. In their model, both historical and current motion vectors are considered as inputs so that the pairwise association between them is captured. However, since the visible layer is in vector form, the spatial information in the multiway data is not retained. In [3] a three-way factored conditional RBM was proposed where a three-way weight tensor is employed to capture the correlations between the input, output, and hidden variables. However, their training data still requires vectorization.

The above works are both aiming to capture the interaction among different vector inputs and are hence not directly applicable to matrix and tensor data. The first RBM designed for tensor inputs is [6], which is called a tensor-variate RBM (TvRBM). In TvRBM, the visible layer is represented as a tensor but the hidden layer is still a vector. Furthermore, the connection between the visible and hidden layers is described by a canonical polyadic (CP) tensor decomposition [9]. However, this CP weight tensor is claimed to constrain the model representation capability [5].

Another RBM related model that utilizes tensor input is the matrix-variate RBM (MvRBM) [5]. The visible and hidden layers in an MvRBM are both matrices. Nonetheless, to limit the number of parameters, an MvRBM models the connection between the visible and hidden layers through two separate matrices, which restricts the ability of the model to capture correlations between different data modes.

All these issues have motivated the MPORBM. Specifically, MPORBM not only employs tensorial visible and hidden layers, but also utilizes a general and powerful tensor network, namely an MPO, to connect the tensorial visible and hidden layers. By doing so, an MPORBM achieves a more powerful model representation capacity than MvRBM and at the same time greatly reduces the model parameters compared to a standard RBM. Note that a mapping of the standard RBM with tensor networks has been described in [10]. However, their work does not generalize the standard RBM to tensorial inputs and is therefore still based on visible and hidden units in vector forms.

Iii Preliminaries

We review some necessary tensor basics, the MPO tensor decomposition, the standard RBM and its tensorial variants.

Iii-a Tensor basics

Tensors are multi-dimensional arrays that are higher-order generalization of vectors (first-order tensors) and matrices (second-order tensors). A th-order (also called -way or -mode) tensor is denoted as where each entry is determined by indices . The numbers are called the dimensions of the tensor. We use boldface capital calligraphic letters , , …to denote tensors, boldface capital letters , , …to denote matrices, boldface letters , , …to denote vectors, and roman letters , , …to denote scalars. and are the transposes of a matrix and a vector . An intuitive and useful graphical representation, named tensor network diagrams, of scalars, vectors, matrices and tensors is depicted in Fig. 1. The unconnected edges are the indices of the tensor. We will mainly employ these graphical representations to visualize the tensor networks and operations in the following sections and refer to [7] for more details. We now briefly introduce some important tensor operations.

Fig. 1: Graphical representation of a scalar , vector , matrix , and third-order tensor .
Fig. 2: Tensor contraction between a third-order tensor and a fourth-order tensor , where the summation runs over and .
Definition 1

(Tensor index contraction): A tensor index contraction is the sum over all possible values of the repeated indices in a set of tensors.

For example, the following contraction between a third-order tensor and a fourth-order tensor

over the and indices results in a third-order tensor . The corresponding tensor network diagram for these index contractions is shown in Fig. 2, where the summation over the and indices is indicated by the connected edges. The resulting diagram has three unconnected indices (), which confirms that is of third order.

Iii-B Tensor MPO decomposition

An MPO is essentially a tensor network representation of a matrix. To relate the row and column indices of a matrix to the corresponding multi-indices of the MPO we need the following definition.

Definition 2

The mapping between a linear index of a vector and its corresponding multi-index when reshaped into a tensor is

Fig. 3: Matrix product operator (MPO) decomposition.

Now suppose we have a matrix , where the index mapping (1) is used to split both the row index and column index into multi-indices , respectively. After this index splitting we can arrange all matrix entries into a -way tensor . Per definition, the corresponding MPO decomposition represents each entry of as


The “building blocks" of the MPO are the -way tensors , also called the MPO-cores. The dimensions of the summation indices are called the MPO-ranks. Observe that the second summation index of is the same as the first index of , which ensures that the right-hand side of (2) results in a scalar. When , the MPO is said to have periodic boundary conditions. We will assume throughout the remainder of this article that . The tensor network diagram representation of (2) is shown in Fig. 3. All index contractions are again represented as edges connecting two tensors. The main motivation to use an MPO representation of a matrix is to reduce its storage requirement. Indeed, using an MPO-representation with maximal MPO-rank for a matrix reduces the storage requirement from down to approximately . This leads to significant storage savings when is relatively small. Consequently, all computations are done on the MPO-cores rather than on the whole matrix itself. This could potentially reduce the computational complexity of the learning algorithm, especially if the visible and hidden layers are also represented as matrix product states (MPS, also called tensor trains in the scientific computing community [11]). This extension is kept for future work.

Iii-C Standard RBM

The standard RBM [1] is a bipartite undirected probabilistic graphical model with one visible layer and one hidden layer , both in vector forms. Here we mainly consider binary RBMs, which implies that entries in both and attain binary values. The standard RBM assigns the following energy function for a specific joint configuration :


where and are the bias of the visible layer and hidden layer, respectively, and is the mapping matrix. All model parameters together are denoted . We can easily use a tensor network diagram to represent Eq. 3 in Fig. 4 (a). The conditional distributions over the hidden and visible layers can be written as:



is the logistic sigmoid function and

denotes a vector of ones. The parameter training in a standard RBM is commonly performed by the contrastive divergence (CD) algorithm [1] and its variants [12, 13].

Iii-D Matrix-variate and tensor-variate RBMs

Here we briefly introduce the two non-vector RBM models, namely MvRBM [5] and TvRBM [6]. TvRBM employs a tensorial visible layer and keeps the hidden layer in a vector form. A rank- CP tensor decomposition is used to connect the visible and hidden layers. However, such a connection form is also criticized to limit the model capability [5]. The corresponding energy function of TvRBM is shown in Fig. 4(b) as a tensor network diagram.

In the MvRBM model, both the input and hidden layers are matrices and are interconnected through two independent weight matrices . This particular construction reduces the total number of parameters from down to but comes at the cost of a limited representation ability, as the weight matrix is constrained to be a Kronecker product of the matrices. The energy function of the MvRBM is graphically represented as a tensor diagram in Fig. 4(c).

Iv Mporbm

We now describe the proposed MPORBM and its customized model parameter learning algorithms.

Iv-a Model definition

In an MPORBM, both the visible layer and the hidden layer are -way tensors. As a result, the weight matrix is now a -way tensor

, which is represented by an MPO instead. With both the visible and hidden layers being tensors, it is therefore also required that the bias vectors

are tensors , respectively. The energy function of an MPORBM is shown in Fig. 4(d) for the specific case where . The vertical edges between the different MPO-cores are the key ingredient in being able to express generic weight tensors , as opposed to the two disconnected weight matrices in the MvRBM model. The corresponding conditional distribution equations over the hidden and visible layers can be derived from the energy function and are visualized in Fig. 4(e)&(f) for the case where . This involves the summation of the weight MPO with either the hidden or visible layer tensors into a -way tensor, which is then added elementwise with the corresponding bias tensor. The final step in the computation of the conditional probability is an elementwise application of the logistic sigmoid function on the resulting tensor.

Fig. 4: Negative energy functions () of (a) RBM; (b) TvRBM; (c) MvRBM; (d) MPORBM. And conditional probability of an MPORBM over its (e) hidden layer; (f) visible layer.

Iv-B MPORBM Learning algorithm

Let denote the model parameters. The model learning task is then formulated into maximizing the training data likelihood:


with respect to model parameter . Similar to the standard RBM [1], the expression of the gradient of the log-likelihood is:


where is the data expectation w.r.t. , which can be computed by Fig. 4(f). is the model expectation w.r.t. , which can be approximately computed by the CD procedure. The main idea in the CD procedure is as follows: first, a Gibbs chain is initialized with one particular training sample = . Figs. 4(e)&(f) are then computed times in an alternating fashion, which results in the chain . The model expectation is then approximated by .

Fig. 5: The derivatives of the log-likelihood function with respect to the th MPO-core .

The derivative of the log-likelihood with respect to the -th MPO-core is visualized in Fig. 5. This involves the computation of 2 fourth-order tensors, obtained by removing the -th MPO-core from the two tensor networks and summing over all connected edges. The resulting tensors are then subtracted elementwise from one another. The derivatives of the log-likelihood with respect to the bias tensors are

We mainly use the CD procedure to train the MPORBM model. However, instead of updating all the MPO-cores simultaneously with one batch of input training data, we employ the alternating optimization procedure (AOP). This involves updating only one MPO-core at a time while keeping the others unchanged using the same batch of input training data. We name this parameter learning algorithm CD-AOP and its pseudo-code is given in Algorithm 1 . The superiority of this alternating optimization procedure over simultaneously updating all MPO-cores, which we call CD-SU from hereon will be demonstrated through numerical experiments.

Algorithm 1

MPORBM learning algorithm (CD-AOP)
Input: Training data of N tensors , the maximum iteration number , the batch size , the momentum , the learning rate and CD step .
Output: Model parameters .

1:  Randomly initialize . Set the bias and the gradient increments .
2:  for iteration number  do
3:     Randomly divide into M batches of size .
4:     for batch  do
5:        for   do
6:           For all data run Gibbs sampling:
7:           for k=0,…,K-1 do
8:              sample according to Fig. 4 (f) with ;
9:              sample according to Fig. 4(e) with
10:           end for
15:        end for
16:     end for
17:  end for

V Experiments

In this section, the MPORBM is compared with both the standard RBM and MvRBM. Specifically, we first investigate the performance and scalability of these models as a means to do feature extraction for classification, together with the influence of the MPO-ranks on the expressive power of the MPORBM. Furthermore, we demonstrate the generative capacity of MPORBM by implementing an image completion task and an image denoising task. For hyperparameter setting, we select mostly the same default values as in MvRBM paper, namely momentum

= , CD step = . For learning rate, we provide a set of possible value and choose the one which achieves the best validation accuracy. Moreover, we choose the same maximum iteration number and batch size for all methods. All experiments are run in MATLAB on an Intel i5 3.2GHz desktop with 16GB RAM.

V-a Data Classification

In the first experiment, we demonstrate the data classification accuracy superiority of MPORBM on extensive standard machine learning datasets, namely the Binary Alphadigits dataset

111 roweis/data.html, normalized DrivFace222, Arcene333 and COIL-100 dataset444 The size of those datasets vary from 320 to 49152. Since we assume binary input in our RBM setting, thus for non-binary datasets, we first use a multi-bit vector to represent each value in the original data. The additional multi-bit dimension is combined with the RGB channel if the dataset is color image. The following table I shows the detailed information of these datasets, while table II and table III show the visual and hidden layer structure of different RBM models for different datasets.

Datasets Original data size Data value range
Alphadigits x Binary
DrivFace x - integer
Arcene x - integer
COIL-100 x x - integer
TABLE I: Detailed information for different datasets.
Alphadigits x x x
DrivFace x x x x
Arcene x x x x
COIL-100 x x x x
TABLE II: Visual layer structure of different RBM models for different datasets.
Alphadigits x x x
DrivFace x x x x
Arcene x x x x
COIL-100 x x x x
TABLE III: Hidden layer structure of different RBM models for different datasets.

We randomly separate the above four datasets into training, validation and testing parts respectively. It might happen that training sample sizes are insufficiently large to train a good RBM model. Thus the training data number we chose for each class is significant smaller than their data dimension, commonly less than 35. The MPORBM-ranks were chosen empirically from a range of 10 to 50 depended on the model input size. The trained RBM models were then employed to extract features from the hidden layer and then those new features will be utilized to train a Nearest Neighbor (

-NN) classifier with

for all experiments. Table IV lists the resulting classification errors.

Alphadigits /
DrivFace /
Arcene /
COIL-100 /
TABLE IV: Classification errors of different RBM models.

The restrictive Kronecker product assumption of the weight matrix in MvRBM explains why it has the worst classification performance in all datasets. The standard RBM also performs worse than MPORBM, which may be caused by overfitting because of the small training sample size. It is notable that due to the PC memory constraint, the standard RBM fails to learn the enormous number of parameters in the full weight matrix for COIL-100 dataset, but MPORBM still gets the best classification performance due to the general MPO weight structure. Moreover, the CD-AOP algorithm shows a higher classification accuracy than CD-SU, which further indicates that the alternating updating scheme is more suitable for our MPORBM model.

V-B Investigation of influence MPO-rank on classification error

In this experiment, we demonstrate how the expressive power of the MPORBM is determined by the MPO-ranks. For this purpose, the MNIST dataset555 was used to train an MPORBM for image classification. The MNIST dataset has a training set of samples and a test set of samples. Each sample is a grayscale picture of handwritten digits {}, where each pixel value was converted into a binary number. The hidden layer states in the RBM model were also regarded as features extracted from the training data and used in a -NN classifier with . The dimensions of the hidden layer were set to . The MPO-rank was varied from 2 up to 200. To reveal the general rule of MPORBM model expression power on MPO-ranks, the above mentioned experiments were run for training sample sizes of 3000, 4000 and 5000, which are randomly chosen from training set. In addition, 10000 samples of the training set were chosen for validation. The classification error for each of these runs is shown in Fig. 6 as a function of . Note that the MPORBM with is identical with the MvRBM for these matrix inputs. It can be seen that the classification error stabilizes when , which indicates that a low-rank MPO may be powerful enough to model the latent probability distribution of the training data and a lot of storage and computational resources can be saved. We need to mention that by setting the MPO-ranks to their maximal values the MPORBM gets the same model expression ability as a standard RBM. This experiments, however, shows that using a low-rank MPORBM is more than enough to model real-life data.

Fig. 6: Classification error as a function of the MPO-rank in the MPORBM model.

V-C Image completion and denoising

In this experiment, we show that an MPORBM is good at generative modeling by implementing image completion and denoising. We test the above generative tasks on the binarized MNIST dataset. Specifically, we randomly choose

samples from the 60k training set to train the RBM models with iterations and using all test samples to check the model generative performance. The standard RBM is set up with vectorial visual layer and hidden layer, while both MvRBM and MPORBM are constructed with 2-way tensorial visual layer and hidden layer. The MPO-rank in this experiment is set to 40 according to the observation from the previous experiment. The visual-hidden layer structure and total weight parameter number of each RBM model is listed in Table V.

Visual structure x x x
Hidden structure x x x
# parameters
TABLE V: The visual-hidden layer structure and total weight parameter number of different RBM models.

In the image completion task, one half of the image is provided to the trained RBM models to complete the another half. As mentioned in [14], RBM completion performance is different when the given half image is in the column or row direction. Thus we investigate the completion ability of different RBM models with giving the right half and the bottom half image respectively. Figure 7 and Figure 8 demonstrate the completed images of different RBM models when given the same randomly selected right and bottom halfs, respectively. It is clear that MvRBM is not able to complete the image, which further confirms the necessity of the MPO generalization. Table VI further lists the average PSNR results between the completed image and the original binarized image on the whole test set. We found that the MPORBM demonstrates a comparable image completion performance to a standard RBM, but with much fewer model parameter (around %).

Fig. 7: Image completion result when given the right half image. First row: original binarized images; second row: RBM completed images; third row: MvRBM completed images; fourth row: MPORBM completed images.
Fig. 8: Image completion result when given the bottom half image. First row: original binarized images; second row: RBM completed images; third row: MvRBM completed images; fourth row: MPORBM completed images.
Right half  dB  dB  dB
Bottom half  dB  dB  dB
TABLE VI: The average PSNR results on the whole test set when given the right and bottom half images respectively, and completing the left and upper half images.

In the image denoising task, we employ the same trained RBM models from the completion task and randomly add a certain percentage of salt & pepper noise to the test set. The average PSNR results between the original binarized pictures and denoised pictures on the whole test set are listed in Table VII. It is clear to see that the denoising performance of the MPORBM is comparable with a standard RBM when % noise is added, but performs better when higher percentages of noise are added. The fewer model parameters in the MPORBM may lead to a more robust generative model.

%  dB  dB  dB
%  dB  dB  dB
%  dB  dB  dB
TABLE VII: The average PSNR results on the whole test set when adding % salt & pepper noises.

Vi Conclusion

This paper has proposed the MPORBM, which preserves the tensorial nature of the input data and utilizes a matrix product operator (MPO) to represent the weight matrix. The MPORBM generalizes all existing RBM models to tensor inputs and has better storage complexity since the number of parameters grows only linearly with the order of the tensor. Furthermore, by representing both the visible and hidden layers as tensors, it is possible to retain useful structural information in the original data. In order to facilitate the exposition on the various computations, tensor network diagrams were introduced and used throughout the paper.

There are several avenues for future research. If both the visible and hidden layers are represented in a matrix product state (MPS) form, then all computations can be done on individual cores. This can significantly improve the computational complexity of the learning algorithm.

Furthermore, analogous to stacking RBMs into a Deep Belief Network (DBN) or Deep Belief Machine (DBM), MPORBMs can be stacked into an MPODBN or MPODBM to extend its expressive power and application. Again, this stacking can also gain computational benefits from keeping all layers in the MPS form.


  • [1] G. E. Hinton, “Training products of experts by minimizing contrastive divergence,” Neural computation, vol. 14, no. 8, pp. 1771–1800, 2002.
  • [2] H. Larochelle and Y. Bengio, “Classification using discriminative restricted Boltzmann machines,” in Proceedings of the 25th international conference on Machine learning.    ACM, 2008, pp. 536–543.
  • [3] A. Krizhevsky, G. Hinton et al., “Factored 3-way restricted Boltzmann machines for modeling natural images,” in

    Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics

    , 2010, pp. 621–628.
  • [4] A.-r. Mohamed and G. Hinton, “Phone recognition using restricted Boltzmann machines,” in Acoustics Speech and Signal Processing (ICASSP), 2010 IEEE International Conference on.    IEEE, 2010, pp. 4354–4357.
  • [5] G. Qi, Y. Sun, J. Gao, Y. Hu, and J. Li, “Matrix variate restricted Boltzmann machine,” in Neural Networks (IJCNN), 2016 International Joint Conference on.    IEEE, 2016, pp. 389–395.
  • [6] T. D. Nguyen, T. Tran, D. Q. Phung, S. Venkatesh et al., “Tensor-Variate Restricted Boltzmann Machines.” in AAAI, 2015, pp. 2887–2893.
  • [7] R. Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Annals of Physics, vol. 349, pp. 117–158, 2014.
  • [8] G. W. Taylor and G. E. Hinton, “Factored conditional restricted Boltzmann machines for modeling motion style,” in Proceedings of the 26th annual international conference on machine learning.    ACM, 2009, pp. 1025–1032.
  • [9] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM review, vol. 51, no. 3, pp. 455–500, 2009.
  • [10] J. Chen, S. Cheng, H. Xie, L. Wang, and T. Xiang, “Equivalence of restricted Boltzmann machines and tensor network states,” Phys. Rev. B, vol. 97, p. 085104, Feb 2018.
  • [11] I. V. Oseledets, “Tensor-train decomposition,” SIAM J. Sci. Comput., vol. 33, no. 5, pp. 2295–2317, 2011.
  • [12] T. Tieleman, “Training restricted Boltzmann machines using approximations to the likelihood gradient,” in Proceedings of the 25th international conference on Machine learning.    ACM, 2008, pp. 1064–1071.
  • [13] T. Tieleman and G. Hinton, “Using fast weights to improve persistent contrastive divergence,” in Proceedings of the 26th Annual International Conference on Machine Learning.    ACM, 2009, pp. 1033–1040.
  • [14] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, “Unsupervised generative modeling using matrix product states,” Physical Review X, vol. 8, no. 3, p. 031012, 2018.