Title: AN INFORMATION-THEORETIC FRAMEWORK FOR FAST AND ROBUST UNSUPERVISED LEARNING VIA NEURAL POPULATION INFOMAX

Abstract: A framework is presented for unsupervised learning of representations based on infomax principle for large-scale neural populations. We use an asymptotic approximation to the Shannon's mutual information for a large neural population to demonstrate that a good initial approximation to the global information-theoretic optimum can be obtained by a hierarchical infomax method. Starting from the initial solution, an efficient algorithm based on gradient descent of the final objective function is proposed to learn representations from the input datasets, and the method works for complete, overcomplete, and undercomplete bases. As confirmed by numerical experiments, our method is robust and highly efficient for extracting salient features from input datasets. Compared with the main existing methods, our algorithm has a distinct advantage in both the training speed and the robustness of unsupervised representation learning. Furthermore, the proposed method is easily extended to the supervised or unsupervised model for training deep structure networks.

1. INTRODUCTION:

How to discover the unknown structures in data is a key task for machine learning. Learning good representations from observed data is important because a clearer description may help reveal the underlying structures. Representation learning has drawn considerable attention in recent years (Bengio et al., 2013) . One category of algorithms for unsupervised learning of representations is based on probabilistic models (Lewicki & Sejnowski, 2000; Hinton & Salakhutdinov, 2006; Lee et al., 2008) , such as maximum likelihood (ML) estimation, maximum a posteriori (MAP) probability estimation, and related methods. Another category of algorithms is based on reconstruction error or generative criterion (Olshausen & Field, 1996; Aharon et al., 2006; Vincent et al., 2010; Mairal et al., 2010; Goodfellow et al., 2014) , and the objective functions usually involve squared errors with additional constraints. Sometimes the reconstruction error or generative criterion may also have a probabilistic interpretation (Olshausen & Field, 1997; Vincent et al., 2010 ). Shannon's information theory is a powerful tool for description of stochastic systems and could be utilized to provide a characterization for good representations (Vincent et al., 2010) . However, computational difficulties associated with Shannon's mutual information (MI) (Shannon, 1948) have hindered its wider applications. The Monte Carlo (MC) sampling (Yarrow et al., 2012) is a convergent method for estimating MI with arbitrary accuracy, but its computational inefficiency makes it unsuitable for difficult optimization problems especially in the cases of high-dimensional input stimuli and large population networks. Bell and Sejnowski (Bell & Sejnowski, 1995; 1997) have directly applied the infomax approach (Linsker, 1988) to independent component analysis (ICA) of data with independent non-Gaussian components assuming additive noise, but their method requires that the number of outputs be equal to the number of inputs. The extensions of ICA to overcomplete or undercomplete bases incur increased algorithm complexity and difficulty in learning of parameters (Lewicki & Sejnowski, 2000; Kreutz-Delgado et al., 2003; Karklin & Simoncelli, 2011) . Since Shannon MI is closely related to ML and MAP (Huang & Zhang, 2017) , the algorithms of representation learning based on probabilistic models should be amenable to information-theoretic treatment. Representation learning based on reconstruction error could be accommodated also by information theory, because the inverse of Fisher information (FI) is the Cramér-Rao lower bound on the mean square decoding error of any unbiased decoder (Rao, 1945) . Hence minimizing the reconstruction error potentially maximizes a lower bound on the MI (Vincent et al., 2010) .

Related problems arise also in neuroscience. It has long been suggested that the real nervous systems might approach an information-theoretic optimum for neural coding and computation (Barlow, 1961; Atick, 1992; Borst & Theunissen, 1999) . However, in the cerebral cortex, the number of neurons is huge, with about 10 5 neurons under a square millimeter of cortical surface (Carlo & Stevens, 2013) . It has often been computationally intractable to precisely characterize information coding and processing in large neural populations.

To address all these issues, we present a framework for unsupervised learning of representations in a large-scale nonlinear feedforward model based on infomax principle with realistic biological constraints such as neuron models with Poisson spikes. First we adopt an objective function based on an asymptotic formula in the large population limit for the MI between the stimuli and the neural population responses (Huang & Zhang, 2017) . Since the objective function is usually nonconvex, choosing a good initial value is very important for its optimization. Starting from an initial value, we use a hierarchical infomax approach to quickly find a tentative global optimal solution for each layer by analytic methods. Finally, a fast convergence learning rule is used for optimizing the final objective function based on the tentative optimal solution. Our algorithm is robust and can learn complete, overcomplete or undercomplete basis vectors quickly from different datasets. Experimental results showed that the convergence rate of our method was significantly faster than other existing methods, often by an order of magnitude. More importantly, the number of output units processed by our method can be very large, much larger than the number of inputs. As far as we know, no existing model can easily deal with this situation.

2. METHODS:



2.1. APPROXIMATION OF MUTUAL INFORMATION FOR NEURAL POPULATIONS:

Suppose the input x is a K-dimensional vector, x = (x 1 , • • • , x K ) T , the outputs of N neurons are denoted by a vector, r = (r 1 , • • • , r N ) T , where we assume N is large, generally N K. We denote random variables by upper case letters, e.g., random variables X and R, in contrast to their vector values x and r. The MI between X and R is defined by I(X; R) = ln p(x|r) p(x) r,x , where

• r,x denotes the expectation with respect to the probability density function (PDF) p(r, x).

Our goal is to maxmize MI I(X; R) by finding the optimal PDF p(r|x) under some constraint conditions, assuming that p(r|x) is characterized by a noise model and activation functions f (x; θ n ) with parameters θ n for the n-th neuron (n = 1, • • • , N ). In other words, we optimize p(r|x) by solving for the optimal parameters θ n . Unfortunately, it is intractable in most cases to solve for the optimal parameters that maximizes I(X; R). However, if p(x) and p(r|x) are twice continuously differentiable for almost every x ∈ R K , then for large N we can use an asymptotic formula to approximate the true value of I(X; R) with high accuracy (Huang & Zhang, 2017) :

EQUATION

where det (•) denotes the matrix determinant and H(X) = -ln p(x) x is the stimulus entropy, G(x) = J(x) + P (x) , (2)

EQUATION

EQUATION

)

Assuming independent noises in neuronal responses, we have p(r|x) = N n=1 p(r n |x; θ n ), and the Fisher information matrix becomes J(x) ≈ N K1 k=1 α k S(x; θ k ), where S(x;

θ k ) = ∂ ln p(r|x;θ k ) ∂x ∂ ln p(r|x;θ k ) ∂x T r|x and α k > 0 (k = 1, • • • , K 1 )

is the population density of parameter θ k , with K1 k=1 α k = 1, and 1 ≤ K 1 ≤ N (see Appendix A.1 for details). Since the cerebral cortex usually forms functional column structures and each column is composed of neurons with the same properties (Hubel & Wiesel, 1962) , the positive integer K 1 can be regarded as the number of distinct classes in the neural population.

Therefore, given the activation function f (x; θ k ), our goal becomes to find the optimal population distribution density α k of parameter vector θ k so that the MI between the stimulus x and the response r is maximized. By Eq. (1), our optimization problem can be stated as follows:

EQUATION

subject to Huang & Zhang, 2017) , we can readily find the optimal solution for small K by efficient numerical methods. For large K, however, finding an optimal solution by numerical methods becomes intractable. In the following we will propose an alternative approach to this problem. Instead of directly solving for the density distribution {α k }, we optimize the parameters {α k } and {θ k } simultaneously under a hierarchical infomax framework.

K1 k=1 α k = 1, α k > 0, ∀k = 1, • • • , K 1 . (6) Since Q G [{α k }] is a convex function of {α k } (

2.2. HIERARCHICAL INFOMAX:

For clarity, we consider neuron model with Poisson spikes although our method is easily applicable to other noise models. The activation function f (x; θ n ) is generally a nonlinear function, such as sigmoid and rectified linear unit (ReLU) (Nair & Hinton, 2010) . We assume that the nonlinear function for the n-th neuron has the following form: f (x; θ n ) = f (y n ; θn ), where

EQUATION

)

with w n being a K-dimensional weights vector, f (y n ; θn ) is a nonlinear function, θ n = (w T n , θT n ) T and θn are the parameter vectors

(n = 1, • • • , N ).

In general, it is very difficult to find the optimal parameters, θ n , n = 1, • • • , N , for the following reasons. First, the number of output neurons N is very large, usually N K. Second, the activation function f (x; θ n ) is a nonlinear function, which usually leads to a nonconvex optimization problem. For nonconvex optimization problems, the selection of initial values often has a great influence on the final optimization results. Our approach meets these challenges by making better use of the large number of neurons and by finding good initial values by a hierarchical infomax method.

We divide the nonlinear transformation into two stages, mapping first from x to y n (n = 1, • • • , N ), and then from y n to f (y n ; θn ), where y n can be regarded as the membrane potential of the n-th neuron, and f (y n ; θn ) as its firing rate. As with the real neurons, we assume that the membrane potential is corrupted by noise:

EQUATION

where Z n ∼ N 0, σ 2 is a normal distribution with mean 0 and variance σ 2 . Then the mean membrane potential of the k-th class subpopulation with N k = N α k neurons is given by

EQUATION

) Define vectors y = (y 1 , • • • , yN ) T , ȳ = (ȳ 1 , • • • , ȳK1 ) T and y = (y 1 , • • • , y K1 ) T , where y k = w T k x (k = 1, • • • , K 1 ). Notice that yn (n = 1, • • • , N

) is also divided into K 1 classes, the same as for r n . If we assume f (x; θ k ) = f (ȳ k ; θk ), i.e. assuming an additive Gaussian noise for y n (see Eq. 9), then the random variables X, Y , Y , Ȳ and R form a Markov chain, denoted by 1 ), and we have the following proposition (see Appendix

X → Y → Y → Ȳ → R (see Figure

EQUATION

EQUATION

and for large

EQUATION

EQUATION

A major advantage of incorporating membrane noise is that it facilitates finding the optimal solution by using the infomax principle. Moreover, the optimal solution obtained this way is more robust; that is, it discourages overfitting and has a strong ability to resist distortion. With vanishing noise

σ 2 → 0, we have Ȳk → Y k , f (ȳ k ; θk ) f (y k ; θk ) = f (x; θ k )

, so that Eqs. ( 13) and ( 14) hold as in the case of large N k .

To optimize MI I(Y ; R), the probability distribution of random variable Y , p(y), needs to be determined, i.e. maximizing I(Y ; R) about p(y) under some constraints should yield an optimal distribution: p * (y) = arg max p(y) I(Y ; R). Let C = max p(y) I (Y ; R) be the channel capacity of neural population coding, and we always have I(X; R) ≤ C (Huang & Zhang, 2017) . To find a suitable linear transformation from X to Y that is compatible with this distribution p * (y), a reasonable choice is to maximize I(X; Y ) (≤ I(X; Y )), where Y is a noise-corrupted version of Y . This implies minimum information loss in the first transformation step. However, there may exist many transformations from X to Y that maximize I(X; Y ) (see Appendix A.3.1). Ideally, if we can find a transformation that maximizes both I(X; Y ) and I(Y ; R) simultaneously, then I(X; R) reaches its maximum value:

I(X; R) = max p(y) I (Y ; R) = C.

From the discussion above we see that maximizing I(X; R) can be divided into two steps, namely, maximizing I(X; Y ) and maximizing I(Y ; R). The optimal solutions of max I(X; Y ) and max I(Y ; R) will provide a good initial approximation that tend to be very close to the optimal solution of max I(X; R).

Similarly, we can extend this method to multilayer neural population networks. For example, a twolayer network with outputs R (1) and R (2) form a Markov chain, 2) , where random variable R( 1) is similar to Y , random variable R (1) is similar to Y , and R( 1) is similar to Ȳ in the above. Then we can show that the optimal solution of max I(X; R (2) ) can be approximated by the solutions of max I(X; R (1) ) and max

X → R(1) → R (1) → R(1) → R (

I( R(1) ; R (2) ), with I( R(1) ; R (2) ) I( R(1) ; R (2) ).

More generally, consider a highly nonlinear feedforward neural network that maps the input x to output z, with z

= F (x;θ) = h L • • • • • h 1 (x), where h l (l = 1, • • • , L

) is a linear or nonlinear function (Montufar et al., 2014) . We aim to find the optimal parameter θ by maximizing I (X; Z). It is usually difficult to solve the optimization problem when there are many local extrema for F (x;θ). However, if each function h l is easy to optimize, then we can use the hierarchical infomax method described above to get a good initial approximation to its global optimization solution, and go from there to find the final optimal solution. This information-theoretic consideration from the neural population coding point of view may help explain why deep structure networks with unsupervised pre-training have a powerful ability for learning representations.

2.3. THE OBJECTIVE FUNCTION:

The optimization processes for maximizing I(X; Y ) and maximizing I(Y ; R) are discussed in detail in Appendix A.3. First, by maximizing I(X; Y ) (see Appendix A.3.1 for details), we can get the optimal weight parameter w k (k = 1, • • • , K 1 , see Eq. 7) and its population density α k (see Eq. 6) which satisfy

EQUATION

EQUATION

where

a = K 1 K -1 0 , C = [c 1 , • • • , c K1 ] ∈ R K0×K1 , CC T = I K0 , I K0 is a K 0 × K 0 identity matrix with integer K 0 ∈ [1, K], the diagonal matrix Σ 0 ∈ R K0×K0 and matrix U 0 ∈ R K×K0 are

given in (A.44) and (A.45), with K 0 given by Eq. (A.52). Matrices Σ 0 and U 0 can be obtained by Σ and U with U T 0 U 0 = I K0 and U 0 Σ 0 U T 0 ≈ UΣU T ≈ xx T x (see Eq. A.23). The optimal weight parameter w k (15) means that the input variable x must first undergo a whiteninglike transformation x = Σ -1/2 0 U T 0 x, and then goes through the transformation y = aC T x, with matrix C to be optimized below. Note that weight matrix W satisfies rank(W) = min(K 0 , K 1 ), which is a low rank matrix, and its low dimensionality helps reduce overfitting during training (see Appendix A.3.1).

By maximizing I (Y ; R) (see Appendix A.3.2), we further solve the the optimal parameters θk for the nonlinear functions f (y k ; θk ), k = 1, • • • , K 1 . Finally, the objective function for our optimization problem (Eqs. 5 and 6) turns into (see Appendix A.3.3 for details):

EQUATION

)

EQUATION

)

where

Φ = diag φ(ŷ 1 ) 2 , • • • , φ(ŷ K1 ) 2 , φ(ŷ k ) = a -1 |∂g k (ŷ k )/∂ ŷk | (k = 1, • • • , K 1 ), g k (ŷ k ) = 2 f (ŷ k ; θk ), ŷk = a -1 y k = c T k x, and x = Σ -1/2 0 U T 0 x.

We apply the gradient descent method to optimize the objective function, with the gradient of Q[C] given by:

EQUATION

where

ω = (ω 1 , • • • , ω K1 ) T , ω k = φ(ŷ k )φ (ŷ k )c T k C ΦC T -1 c k , k = 1, • • • , K 1 .

When K 0 = K 1 (or K 0 > K 1 ), the objective function Q[C] can be reduced to a simpler form, and its gradient is also easy to compute (see Appendix A.4.1). However, when K 0 < K 1 , it is computationally expensive to update C by applying the gradient of Q[C] directly, since it requires matrix inversion for every x. We use another objective function Q[C] (see Eq. A.118) which is an approximation to Q[C], but its gradient is easier to compute (see Appendix A.4.2). The function

Q[C] is the approximation of Q[C],

ideally they have the same optimal solution for the parameter C.

Usually, for optimizing the objective in Eq. 17, the orthogonality constraint (Eq. 18) is unnecessary. However, this orthogonality constraint can accelerate the convergence rate if we employ it for the initial iteration to update C (see Appendix A.5).

3. EXPERIMENTAL RESULTS:

We have applied our methods to the natural images from Olshausen's image dataset (Olshausen & Field, 1996) T . We obtained the initial values (see Appendix A.3.2): b 0 = 0 and β 0 ≈ 1.81 K 1 K -1 0 . For our experiments, we set β = 0.5β 0 for iteration epoch t = 1, • • • , t 0 and β = β 0 for t = t 0 + 1, • • • , t max , where t 0 = 50.

Firstly, we tested the case of K = K 0 = K 1 = 144 and randomly sampled M = 10 5 image patches with size 12×12 from the Olshausen's natural images, assuming that N = 10 6 neurons were divided into K 1 = 144 classes and = 1 (see Eq. A.52 in Appendix). The input patches were preprocessed by the ZCA whitening filters (see Eq. A.68). To test our algorithms, we chose the batch size to be equal to the number of training samples M , although we could also choose a smaller batch size. We updated the matrix C from a random start, and set parameters t max = 300, v 1 = 0.4, and τ = 0.8 for all experiments.

In this case, the optimal solution C looked similar to the optimal solution of IICA (Bell & Sejnowski, 1997) . We also compared with the fast ICA algorithm (FICA) (Hyvärinen, 1999) , which is faster than IICA. We also tested the restricted Boltzmann machine (RBM) (Hinton et al., 2006) for a unsupervised learning of representations, and found that it could not easily learn Gabor-like filters from Olshausen's image dataset as trained by contrastive divergence. However, an improved method by adding a sparsity constraint on the output units, e.g., sparse RBM (SRBM) (Lee et al., 2008) or sparse autoencoder (Hinton, 2010) , could attain Gabor-like filters from this dataset. Similar results with Gabor-like filters were also reproduced by the denoising autoencoders (Vincent et al., 2010) , which method requires a careful choice of parameters, such as noise level, learning rate, and batch size.

In order to compare our methods, i.e. Algorithm 1 (Alg.1, see Appendix A.4.1) and Algorithm 2 (Alg.2, see Appendix A.4.2), with other methods, i.e. IICA, FICA and SRBM, we implemented these algorithms using the same initial weights and the same training data set (i.e. 10 5 image patches preprocessed by the ZCA whitening filters). To get a good result by IICA, we must carefully select the parameters; we set the batch size as 50, the initial learning rate as 0.01, and final learning rate as 0.0001, with an exponential decay with the epoch of iterations. IICA tends to have a faster convergence rate for a bigger batch size but it may become harder to escape local minima. For FICA, we chose the nonlinearity function f (u) = log cosh(u) as contrast function, and for SRBM, we set the sparseness control constant p as 0.01 and 0.03. The number of epoches for iterations was set to 300 for all algorithms. Figure 2 shows the filters learned by our methods and other methods. Each filter in Figure 2 (a) corresponds to a column vector of matrix Č (see Eq. A.69), where each vector for display is normalized by čk Figure 3 shows how the coefficient entropy (CFE) (see Eq. A.122) and the conditional entropy (CDE) (see Eq. A.125) varied with training time. We calculated CFE and CDE by sampling once every 10 epoches from a total of 300 epoches. These results show that our algorithms had a fast convergence rate towards stable solutions while having CFE and CDE values similar to the algorithm of IICA, which converged much more slowly. Here the values of CFE and CDE should be as small as possible for a good representation learned from the same data set. Here we set epoch number t 0 = 50 in our algorithms (see Alg.1 and Alg.2), and the start time was set to 1 second. This explains the step seen in Figure 3 (b) for Alg.1 and Alg.2 since the parameter β was updated when epoch number t = t 0 . FICA had a convergence rate close to our algorithms but had a big CFE, which is reflected by the quality of the filter results in Figure 2 . The convergence rate and CFE for SRBM were close to IICA, but SRBM had a much bigger CDE than IICA, which implies that the information had a greater loss when passing through the system optimized by SRBM than by IICA or our methods.

← čk / max(|č 1,k |, • • • , |č K,k |), k = 1, • • • , K 1 . The results in Figures 2(a

From Figure 3 (c) we see that the CDE (or MI I(X; R), see Eq. A.124 and A.125) decreases (or increases) with the increase of the value of the sparseness control constant p. Note that a smaller p means sparser outputs. Hence, in this sense, increasing sparsity may result in sacrificing some information. On the other hand, a weak sparsity constraint may lead to failure of learning Gaborlike filters (see Figure 2 (f)), and increasing sparsity has an advantage in reducing the impact of noise in many practical cases. Similar situation also occurs in sparse coding (Olshausen & Field, 1997) , which provides a class of algorithms for learning overcomplete dictionary representations of the input signals. However, its training is time consuming due to its expensive computational cost, although many new training algorithms have emerged (e.g. Aharon et al., 2006; Elad & Aharon, 2006; Lee et al., 2006; Mairal et al., 2010) . See Appendix A.5 for additional experimental results.

4. CONCLUSIONS:

In this paper, we have presented a framework for unsupervised learning of representations via information maximization for neural populations. Information theory is a powerful tool for machine learning and it also provides a benchmark of optimization principle for neural information processing in nervous systems. Our framework is based on an asymptotic approximation to MI for a large-scale neural population. To optimize the infomax objective, we first use hierarchical infomax to obtain a good approximation to the global optimal solution. Analytical solutions of the hierarchical infomax are further improved by a fast convergence algorithm based on gradient descent. This method allows us to optimize highly nonlinear neural networks via hierarchical optimization using infomax principle.

From the viewpoint of information theory, the unsupervised pre-training for deep learning (Hinton & Salakhutdinov, 2006; Bengio et al., 2007) may be reinterpreted as a process of hierarchical infomax, which might help explain why unsupervised pre-training helps deep learning (Erhan et al., 2010) . In our framework, a pre-whitening step can emerge naturally by the hierarchical infomax, which might also explain why a pre-whitening step is useful for training in many learning algorithms (Coates et al., 2011; Bengio, 2012) .

Our model naturally incorporates a considerable degree of biological realism. It allows the optimization of a large-scale neural population with noisy spiking neurons while taking into account of multiple biological constraints, such as membrane noise, limited energy, and bounded connection weights. We employ a technique to attain a low-rank weight matrix for optimization, so as to reduce the influence of noise and discourage overfitting during training. In our model, many parameters are optimized, including the population density of parameters, filter weight vectors, and parameters for nonlinear tuning functions. Optimizing all these model parameters could not be easily done by many other methods.

Our experimental results suggest that our method for unsupervised learning of representations has obvious advantages in its training speed and robustness over the main existing methods. Our model has a nonlinear feedforward structure and is convenient for fast learning and inference. This simple and flexible framework for unsupervised learning of presentations should be readily extended to training deep structure networks. In future work, it would interesting to use our method to train deep structure networks with either unsupervised or supervised learning.

Proposition 2. The following inequations hold,

1 2 ln det C ΦC T x ≤ 1 2 ln det C Φ x C T , (A.97) ln det CΦC T x ≤ ln det C Φ x C T (A.98) ≤ 1 2 ln det C Φ 2 x C T (A.99) ≤ 1 2 ln det C Φ x C T , (A.100) ln det CΦC T ≤ 1 2 ln det C ΦC T , (A.101)

where C ∈ R K0×K1 , K 0 ≤ K 1 , and CC T = I K0 .

Proof. Functions ln det C Φ x C T and ln det C Φ x C T are concave functions about p (x) (see the proof of Proposition 5.2. in Huang & Zhang, 2017) , which fact establishes inequalities (A.97) and (A.98).

Next we will prove the inequality (A.101). By SVD, we have

CΦ = Ü D VT , (A.102)

where

Ü is a K 0 × K 0 unitary orthogonal matrix, V = [v 1 , v2 , • • • , vK1

] is an K 1 × K 1 unitary orthogonal matrix, and D is an K 0 × K 1 rectangular diagonal matrix with K 0 positive real numbers on the diagonal. By the matrix Hadamard's inequality and Cauchy-Schwarz inequality we have

det CΦC T CΦC T det C ΦC T -1 = det D VT C T C V DT D DT -1 = det VT 1 C T C V1 = det C V1 2 ≤ K0 k=1 C V1 2 k,k ≤ K0 k=1 CC T 2 k,k VT 1 V1 2 k,k = 1, (A.103) where V1 = [v 1 , v2 , • • • , vK0 ] ∈ R K1×K0 .

The last equality holds because of CC T = I K0 and VT 1 V1 = I K0 . This establishes inequality (A.101) and the equality holds if and only if

K 0 = K 1 or C V1 = I K0 .

Similarly, we get inequality (A.99):

ln det C Φ x C T ≤ 1 2 ln det C Φ 2 x C T . (A.104)

By Jensen's inequality, we have

φ (ŷ k ) 2 x ≤ φ (ŷ k ) 2 x , ∀k = 1, • • • , K 1 . (A.105)

Then it follows from (A.105) that inequality (A.100) holds: .106) where x is defined by Eq. (A.68), and wk is the corresponding optimal filter. To estimate the probability density of coefficients q k (y k ) (k = 1, • • • , K 1 ) from the M training samples, we apply the kernel density estimation for q k (y k ) and use a normal kernel with an adaptive optimal window width. Then we define the CFE h as A.123) where q k (y k ) is quantized as discrete q k (n∆) and ∆ is the step size.

EQUATION

EQUATION

Methods such as IICA and SRBM as well as our methods have feedforward structures in which information is transferred directly through a nonlinear function, e.g., the sigmoid function. We can use the amount of transmitted information to measure the results learned by these methods. Consider a neural population with N neurons, which is a stochastic system with nonlinear transfer functions. We chose a sigmoidal transfer function and Gaussian noise with standard deviation set to 1 as the system noise. In this case, from ( 1 

I(X; R) = I X; R = H( X) -H X|R ĨG = H( X) -h 1 , (A.124) H X|R h 1 = - 1 2 ln det 1 2πe N K -1 0 C ΦC T + I K0 x , (A.125)

where we set N = 10 6 . A good representation should make the MI I(X; R) as big as possible.

Equivalently, for the same inputs, a good representation should make the conditional entropy (CDE) H X|R (or h 1 ) as small as possible. A.5.4 IMAGE DENOISING Similar to the sparse coding method applied to image denoising (Elad & Aharon, 2006) , our method (see Eq. A.130) can also be applied to image denoising, as shown by an example in Figure 7 . The filters or bases were learned by using 7×7 image patches sampled from the left half of the image, and subsequently used to reconstruct the right half of the image which was distorted by Gaussian noise.

A common practice for evaluating the results of image denoising is by looking at the difference between the reconstruction and the original image. If the reconstruction is perfect the difference should look like Gaussian noise. In Figure 7 (c) and 7(d) a dictionary (100 bases) was learned by MBDL and orthogonal matching pursuit was used to estimate the sparse solution.1 For our method (shown in Figure 7 (b)), we first get the optimal filters parameter W, a low rank matrix (K 0 < K), then from the distorted image patches x m (m = 1, • • • , M ) we get filter outputs y m = W T x m and the reconstruction xm = By m (parameters: = 0.975 and K 0 = K 1 = 14). As can be seen from Figure 7 , our method worked better than dictionary learning, although we only used 14 bases compared with 100 bases used by dictionary learning. Our method is also more efficient. We can get better optimal bases B by a generative model using our infomax approach (details not shown). 

:

Python source code is available at http://scikit-learn.org/stable/ downloads/plot image denoising.py

Figures and Tables:

and the images of handwritten digits from MNIST dataset(LeCun et al., 1998) using Matlab 2016a on a computer with 12 Intel CPU cores (2.4 GHz). The gray level of each raw image was normalized to the range of 0 to 1. M image patches with size w × w = K for training were randomly sampled from the images. We used the Poisson neuron model with a modified sigmoidal tuning function f (y; θ) = 1 4(1+exp(-βy-b)) 2 , with g(y) = 2 f (y; θ) = 1 1+exp(-βy-b) , where θ = (β, b)

), 2(b) and 2(c) look very similar to one another, and slightly different from the results in Figure2(d) and 2(e). There are no Gabor-like filters in Figure2(f), which corresponds to SRBM with p = 0.03.

Figure 2: Comparison of filters obtained from 10 5 natural image patches of size 12×12 by our methods (Alg.1 and Alg.2) and other methods. The number of output filters was K 1 = 144. (a): Alg.1. (b): Alg.2. (c): IICA. (d): FICA. (e): SRBM (p = 0.01). (f): SRBM (p = 0.03).

Figure 4: Comparison of basis vectors obtained by our method and other methods. Panel (a)-(e) correspond to panel (a)-(e) in Figure 2, where the basis vectors are given by (A.130). The basis vectors in panel (f) are learned by MBDL and given by (A.127).

Figure 5: Filters and bases obtained from Olshausen's image dataset and MNIST dataset by Algorithm 2. (a) and (b): 400 typical filters and the corresponding bases obtained from Olshausen's image dataset, where K 0 = 82 and K 1 = 1024. (c) and (d): 400 typical filters and the corresponding bases obtained from the MNIST dataset, where K 0 = 183 and K 1 = 1024.

Figure 6 shows that CFE as a function of training time for Alg.2, where Figure 6(a) corresponds to Figure 5(a)-5(b) for learning nature image patches and Figure 6(b) corresponds to Figure 5(c)-5(d) for learning MNIST dataset. We set parameters t max = 100 and τ = 0.8 for all experiments and varied parameter v 1 for each experiment, with v 1 = 0.2, 0.4, 0.6 or 0.8. These results indicate a fast convergence rate for training on different datasets. Generally, the convergence is insensitive to the change of parameter v 1 .

Figure 7: Image denoising. (a): the right half of the original image is distorted by Gaussian noise and the norm of the difference between the distorted image and the original image is 23.48. (b): image denoising by our method (Algorithm 1), with 14 bases used. (c) and (d): image denoising using dictionary learning, with 100 bases used.

), (A.8) and (A.11), we see that the approximate MI I G is equivalent to the case of the Poisson neuron model. It follows from (A.70)-(A.82) that

Table content: <table/>

Peer review (iclr2017):
This is an 18 page paper plus appendix which presents a mathematical derivation for infomax for an actual neural population with noise. The original Bell & Sejnowski infomax framework only considered the no noise case. Results are shown for natural image patches and the mnist dataset, which qualitatively resemble results obtained with other methods. This seems like an interesting and potentially more general approach to unsupervised learning. However the paper is quite long and it was difficult for me to follow all the twists and turns. For example the introduction of the hierarchical model was confusing and it took several iterations to understand where this was going. *Hierarchical* is probably not the right terminology here because it*s not like a deep net hierarchy, it*s just decomposing the tuning curve function into different parts. I would recommend that the authors try to condense the paper so that the central message and important steps are conveyed in short order, and then put the more complete mathematical development into a supplementary document. Also, the authors should look at the work of Karklin & Simoncelli 2011 which is highly related. They also use an infomax framework for a noisy neural population to derive on and off cells in the retina, and they show the conditions under which orientation selectivity emerges.