Sparse Activity and Sparse Connectivity in Supervised Learning

by   Markus Thom, et al.
Universität Ulm

Sparseness is a useful regularizer for learning in a wide range of applications, in particular in neural networks. This paper proposes a model targeted at classification tasks, where sparse activity and sparse connectivity are used to enhance classification capabilities. The tool for achieving this is a sparseness-enforcing projection operator which finds the closest vector with a pre-defined sparseness for any given vector. In the theoretical part of this paper, a comprehensive theory for such a projection is developed. In conclusion, it is shown that the projection is differentiable almost everywhere and can thus be implemented as a smooth neuronal transfer function. The entire model can hence be tuned end-to-end using gradient-based methods. Experiments on the MNIST database of handwritten digits show that classification performance can be boosted by sparse activity or sparse connectivity. With a combination of both, performance can be significantly better compared to classical non-sparse approaches.


Finding trainable sparse networks through Neural Tangent Transfer

Deep neural networks have dramatically transformed machine learning, but...

Shoestring: Graph-Based Semi-Supervised Learning with Severely Limited Labeled Data

Graph-based semi-supervised learning has been shown to be one of the mos...

On improving deep learning generalization with adaptive sparse connectivity

Large neural networks are very successful in various tasks. However, wit...

Non-Differentiable Supervised Learning with Evolution Strategies and Hybrid Methods

In this work we show that Evolution Strategies (ES) are a viable method ...

Provable Methods for Training Neural Networks with Sparse Connectivity

We provide novel guaranteed approaches for training feedforward neural n...

Enhancing Explainability of Neural Networks through Architecture Constraints

Prediction accuracy and model explainability are the two most important ...

Efficient combination of pairswise feature networks

This paper presents a novel method for the reconstruction of a neural ne...

1 Introduction

Sparseness is a concept of efficiency in neural networks, and exists in two variants in that context (Laughlin and Sejnowski, 2003). The sparse activity property means that only a small fraction of neurons is active at any time. The sparse connectivity property means that each neuron is connected to only a limited number of other neurons. Both properties have been observed in mammalian brains (Hubel and Wiesel, 1959; Olshausen and Field, 2004; Mason et al., 1991; Markram et al., 1997)

and have inspired a variety of machine learning algorithms. A notable result was achieved through the sparse coding model of

Olshausen and Field (1996). Given small patches from images of natural scenes, the model is able to produce Gabor-like filters, resembling properties of simple cells found in mammalian primary visual cortex (Hubel and Wiesel, 1959; Vinje and Gallant, 2000). Another example is the optimal brain damage method of LeCun et al. (1990), which can be used to prune synaptic connections in a neural network, making connectivity sparse. Although only a small fraction of possible connections remains after pruning, this is sufficient to achieve equivalent classification results. Since then, numerous approaches on how to measure sparseness have been proposed, see Hurley and Rickard (2009) for an overview, and how to achieve sparse solutions of classical machine learning problems.

The pseudo-norm is a natural sparseness measure. Its computation consists of counting the number of non-vanishing entries in a vector. Using it rather than other sparseness measures has been shown to induce biologically more plausible properties (Rehn and Sommer, 2007). However, finding of optimal solutions subject to the pseudo-norm turns out to be NP-hard (Natarajan, 1995; Weston et al., 2003). Analytical properties of this counting measure are very poor, for it is non-continuous, rendering the localization of approximate solutions difficult. The Manhattan norm of a vector is a convex relaxation of the pseudo-norm (Donoho, 2006), and has been employed in a vast range of applications. This sparseness measure has the significant disadvantage of not being scale-invariant, so that an intuitive notion of sparseness cannot be derived from it.

1.1 Hoyer’s Normalized Sparseness Measure

A normalized sparseness measure based on the ratio of the or Manhattan norm and the or Euclidean norm of a vector has been proposed by Hoyer (2004),

where higher values indicate more sparse vectors. is well-defined because holds for all (Laub, 2004). As for all and all , is also scale-invariant. As composition of differentiable functions, is differentiable on its entire domain.

This sparseness measure fulfills all criteria of Hurley and Rickard (2009) except for Dalton’s fourth law, which states that the sparseness of a vector should be identical to the sparseness of the vector resulting from multiple concatenation of the original vector. This property, however, is not crucial for a proper sparseness measure. For example, sparseness of connectivity in a biological brain increases quickly with its volume, so that connectivity in a human brain is about  times more sparse than in a rat brain (Karbowski, 2003). It follows that features all desirable properties of a proper sparseness measure.

A sparseness-enforcing projection operator, suitable for projected gradient descent algorithms, was proposed by Hoyer (2004) for optimization with respect to . For a pre-defined target degree of sparseness , the operator finds the closest vector of a given scale that has sparseness given an arbitrary vector. This can be expressed formally as Euclidean projection onto parameterizations of the sets

The first set is for achieving unrestricted projections, whereas the latter set is useful in situations where only non-negative solutions are feasible, for example in non-negative matrix factorization problems. The constants are target norms and can be chosen such that all points in these sets achieve a sparseness of . For example, if was set to unity for yielding normalized projections, then can be easily derived from the definition of .

Hoyer’s original algorithm for computation of such a projection is an alternating projection onto a hyperplane representing the

norm constraint, a hypersphere representing the norm constraint, and the non-negative orthant. A slightly modified version of this algorithm has been proved to be correct by Theis et al. (2005) in the special case when exactly one negative entry emerges that is zeroed out in the orthant projection. However, there is still no mathematically satisfactory proof for the general case.

1.2 Contributions of this Paper

This paper improves upon previous work in the following ways. Section 2 proposes a simple algorithm for carrying out sparseness-enforcing projections with respect to Hoyer’s sparseness measure. Further, an improved algorithm is proposed and compared with Hoyer’s original algorithm. Because the projection itself is differentiable, it is the ideal tool for achieving sparseness in gradient-based learning. This is exploited in Section 3

, where the sparseness projection is used to obtain a classifier that features both sparse activity and sparse connectivity in a natural way. The benefit of these two key properties is demonstrated on a real-world classification problem, proving that sparseness acts as regularizer and improves classification results. The final sections give an overview of related concepts and conclude this paper.

On the theoretical side, a first rigorous and mathematically satisfactory analysis of the properties of the sparseness-enforcing projection is provided. This is lengthy and technical and therefore deferred into several appendixes. Appendix A fixes the notation and gives an introduction to general projections. In Appendix B, certain symmetries of subsets of the Euclidean space and their effect on projections onto such sets is studied. The problem of finding projections onto sets where Hoyer’s sparseness measure attains a constant value is addressed in Appendix C. Ultimately, the algorithms proposed in Section 2 are proved to be correct. Appendix D investigates analytical properties of the sparseness projection and concludes with an efficient algorithm that computes its gradient. The gradients for optimization of the parameters of the architecture proposed in Section 3 are collected in the final Appendix E.

2 Algorithms for the Sparseness-Enforcing Projection Operator

The projection onto a set is a fundamental concept, for example see Deutsch (2001): Let and . Then every point in

is called Euclidean projection of onto . When there is exactly one point in , then is used as an abbreviation. Because is finite-dimensional, is nonempty for all if and only if is closed, and is a singleton for all if and only if is closed and convex (Deutsch, 2001). In the literature, the elements from are also called best approximations to from .

Projections onto sets that fulfill certain symmetries are of special interest in this paper and are formalized and discussed in Appendix B in greater detail. It is notable that projections onto a permutation-invariant set , that is a set where membership is stable upon coordinate permutation, are order-preserving. This is proved in Lemma B1. As a consequence, when a vector is sorted in ascending or descending order, then its projection onto is sorted accordingly. If is reflection-invariant, that is when the signs of arbitrary coordinates can be swapped without violating membership in , then the projection onto is orthant-preserving, as shown in Lemma B2. This means that a point and its projection onto are located in the same orthant. By exploiting this property, projections onto can be yielded by recording and discarding the signs of the coordinates of the argument, projecting onto , and finally restoring the signs of the coordinates of the result using the signs of the argument. This is formalized in Lemma B.

As an example for these concepts, consider the set of all vectors with exactly non-vanishing entries. is clearly both permutation-invariant and reflection-invariant. Therefore, the projection with respect to an pseudo-norm constraint must be both order-preserving and orthant-preserving. In fact, the projection onto consists simply of zeroing out all entries but the that are greatest in absolute value (Blumensath and Davies, 2009). This trivially fulfills the aforementioned properties of order-preservation and orthant-preservation.

Permutation-invariance and reflection-invariance are closed under intersection and union operations. Therefore, the unrestricted target set for the projection is permutation-invariant and reflection-invariant. It is hence enough to handle projections onto in the first place, as projections onto the unrestricted target set can easily be recovered.

In the remainder of this section, let be the problem dimensionality and let be the fixed target norms, which must fulfill to avoid the existence of only trivial solutions. In the applications of the sparseness projection in this paper, is always set to unity to achieve normalized projections, and is adjusted as explained in Section 1.1 to achieve the target degree of sparseness . The related problem of finding the best approximation to a point regardless of the concrete scaling, that is computing projections onto , can be solved by projecting onto and rescaling the result such as to minimize under variation of , which yields . This method is justified theoretically by Remark 2.

2.1 Alternating Projections

First note that the target set can be written as an intersection of simpler sets. Let be the canonical basis of the -dimensional Euclidean space . Further, let be the vector where all entries are identical to unity. Then denotes the target hyperplane where the coordinates of all points sum up to . In the non-negative orthant , this is equivalent to the norm constraint. Further, define as the target hypersphere of all points satisfying the norm constraint. This yields the following factorization:

For computation of projections onto an intersection of a finite number of closed and convex sets, it is enough to perform alternating projections onto the members of the intersection (Deutsch, 2001). As is clearly non-convex, this general approach has to be altered to work in this specific setup.

First, consider , which denotes the intersection of the norm target hyperplane and the norm target hypersphere. essentially possesses the structure of a hypercircle, that is, all points in lie also in and there is a central point and a real number such that all points in have squared distance from . It will be shown in Appendix C that and . The intersection of the non-negative orthant with the norm hyperplane, , is a scaled canonical simplex. Its barycenter coincides with the barycenter of . Finally, for an index set let denote the subset of points from , where all coordinates with index not in vanish. Its barycenter is given by . With these preparations, a simple algorithm can be proposed; it computes the sparseness-enforcing projection with respect to a constraint induced by Hoyer’s sparseness measure .

Input: and with .
Output: where .
// Project onto target hyperplane and target hypercircle .
1 ;
2 ;
// Perform alternating projections until feasible solution is found.
3 while do
      // Project onto scaled canonical simplex .
4       ;
       // Project onto keeping already vanished coordinates at zero.
5       where ;
Algorithm 1 Proposed algorithm for computing the sparseness-enforcing projection operator for Hoyer’s sparseness measure .

For every , Algorithm 1 computes an element from . If after line 1 and after line 1 in all iterations, then is a singleton. As already pointed out, the idea of Algorithm 1 is that projections onto can be computed by alternating projections onto the geometric structures just defined. The rigorous proof of correctness from Appendix C proceeds by showing that the set of solutions is not tampered by projection onto the intermediate structures , , and . Because of the non-convexity of and , the relation between these sets and the simplex is non-trivial and needs long arguments to be described further, see especially Lemma C.2.2 and Corollary 10.

The projection onto the hyperplane is straightforward and discussed in Section C.1.1. As is essentially a hypersphere embedded in a subspace of , projections of points from onto are achieved by shifting and scaling, see Section C.1.2. The alternating projection onto and in the beginning of Algorithm 1 make the result of the projection onto invariant to positive scaling and arbitrary shifting of the argument, as shown in Corollary C.1.2. This is especially useful in practice, alleviating the need for certain pre-processing methods. The formula for projections onto can be generalized for projections onto for an index set , by keeping already vanished coordinates at zero, see Section C.3.

Projections onto the simplex are more involved and discussed at length in Section C.2. The most relevant result is that if , then there exists a separator such that , where the maximum is taken element-wise (Chen and Ye, 2011). In the cases considered in this paper it is always as shown in Lemma C.3. This implies that all entries in that are less than do not survive the projection, and hence the pseudo-norm of is strictly greater than that of . The simplex projection therefore enhances sparseness.

Input: and .
Output: such that and .
// Sort the input vector in descending order.
1 Let such that and ;
// Find the only feasible separator .
2 ;
3 for to do
4      ; ;
5       if then return ;
7end ; ; return ;
Algorithm 2 Computation of information for performing projections onto , which is a scaled canonical simplex. This is an adapted version of the algorithm of Chen and Ye (2011).
Input: and with .
Output: where .
1 procedure proj_L
        // Compute squared radius of (Lemma C.1.2).
        // (Remark C.1.1).
2      if then
              // equals the barycenter of ,
              // pick a sorted projection (Remark C.1.2).
      else ;
        // Pick unique projection (Lemma C.1.2).
// Beginning of main body.
Let such that and ;
  // Sort the input vector.
  // Project onto (Lemma C.1.1).
  // Project in-place onto .
// Perform alternating projections until feasible solution is found.
  // Store current number of relevant entries of .
6while do
        // This is carried out by Algorithm 2.
        // Project onto (Proposition C.2.1).
        // Project onto where (Lemma C.3).
end // Undo sorting permutation and set remaining entries to zero.
8 ; for to do ;
Algorithm 3 Explicit and optimized variant of Algorithm 1.

The separator and the number of nonzero entries in the projection onto can be computed with Algorithm 2, which is an adapted version of the algorithm of Chen and Ye (2011). In line 2, denotes the symmetric group and denotes the permutation matrix associated with a permutation . The algorithm works by sorting its argument and then determining as the mean value of the largest entries of minus the target norm . The number of relevant entries for computation of is equal to the pseudo-norm of the projection and is found by trying all feasible values, starting with the largest ones. The computational complexity of Algorithm 2 is dominated by sorting the input vector and is thus quasilinear.

2.2 Optimized Variant

Because of the permutation-invariance of the sets involved in the projections, it is enough to sort the vector that is to be projected onto once. This guarantees that the working vector that emerges from subsequent projections is sorted also. No additional sorting has then to be carried out when using Algorithm 2 for projections onto

. This additionally has the side effect that the non-vanishing entries of the working vector are always concentrated in its first entries. Hence all relevant information can always be stored in a small unit-stride array, to which access is more efficient than to a large sparse array. Further, the index set

of non-vanishing entries in the working vector is always of the form , where is the number of nonzero entries.

Algorithm 3 is a variant of Algorithm 1 where these optimizations were applied, and where the explicit formulas for the intermediate projections were used. The following result, which is proved in Appendix C, states that both algorithms always compute the same result: Algorithm 1 is equivalent to Algorithm 3. Projections onto increase the amount of vanishing entries in the working vector, which is of finite dimension . Hence, at most alternating projections are carried out, and the algorithm terminates in finite time. Further, the complexity of each iteration is at most linear in the pseudo-norm of the working vector. The theoretic overall computational complexity is thus at most quadratic in problem dimensionality .

2.3 Comparison with Hoyer’s Original Algorithm

The original algorithm for the sparseness-enforcing projection operator proposed by Hoyer (2004) is hard to understand, and correctness has been proved by Theis et al. (2005) in a special case only. A simple alternative has been proposed with Algorithm 1 in this paper. Based on the symmetries induced by Hoyer’s sparseness measure and by exploiting the projection onto a simplex, an improved method was given in Algorithm 3.

The improved algorithm proposed in this paper always requires at most the same number of iterations of alternating projections as the original algorithm. The original algorithm uses a projection onto the non-negative orthant to achieve vanishing coordinates in the working vector. This operation can be written as . In the improved algorithm, a simplex projection is used for this purpose, expressed formally as with chosen accordingly. Due to the theoretical results on simplex geometry from Section C.2 and their application in Lemma C.3 in Section C.3, the number is always non-negative. Therefore, at least the same amount of entries is set to zero in the simplex projection compared to the projection onto the non-negative orthant, see also Corollary C.3. Hence with induction for the number of non-vanishing entries in the working vector, the number of iterations the proposed algorithm needs to terminate is bounded by the number of iterations the original method needs to terminate given the same input.

The experimental determination of an estimate of the number of iterations required was carried out as follows. Random vectors with sparseness

were sampled and their sparse projections were computed using the respective algorithms, to gain the best normalized approximations with a target sparseness degree of . For both algorithms the very same vectors were used as input. During the run-time of the algorithms, the number of iterations that were necessary to compute the result were counted. Additionally, the number of nonzero entries in the working vector was recorded in each iteration. This was done for different dimensionalities, and for each dimensionality  vectors were sampled.

Figure 1: Comparison of the number of iterations of the original algorithm for the projection onto with the improved version as proposed in this paper. The sparseness-enforcing projection with target sparseness was carried out for input vectors of sparseness . The thick lines indicate the mean number of iterations required for the projection, and the thin lines indicate the minimum and maximum number of iterations, respectively. Even for input vectors with a million entries, less than  iterations are required to find the projection. With the improved algorithm, this reduces to at most  iterations.
Figure 2: Comparison of the number of non-vanishing entries in the working vectors of the original algorithm and the improved algorithm during run-time. The algorithms were run with input vectors of dimensionality and initial sparseness to compute projections with sparseness

. Standard deviations were always less than

; they were omitted in the plot to avoid clutter. The algorithm proposed in this paper reduces dimensionality more quickly and terminates earlier than the original algorithm.
Figure 3: Ratio of the computation time of the original algorithm and the improved algorithm for a variety of input dimensionality and initial vector sparseness. Numbers greater than one indicate parameterizations where the proposed algorithm is more efficient than the original one. There is a large region where the speedup is decent.

Figure 1 shows statistics on the number of iterations the algorithms needed to terminate. As was already observed by Hoyer (2004), the number of required iterations grows very slowly with problem dimensionality. For , only between  and  iterations were needed with the original algorithm to compute the result. With Algorithm 3, this can be improved to requiring  to  iterations, which amounts to roughly  less iterations. Due to the small slope in the number of required iterations, it can be conjectured that this quantity is at most logarithmic in problem dimensionality . If this applies, the complexity of Algorithm 3 is at most quasilinear. Because the input vector is sorted in the beginning, it is also not possible to fall below this complexity class.

The progress of working dimensionality reduction for problem dimensionality is depicted in Figure 2, averaged over the  input vectors from the experiment. After the first iteration, that is after projecting onto and , the working dimensionality still matches the input dimensionality. Starting with the second iteration, dimensions are discarded by projecting onto in the original algorithm and onto in the improved variant, which yields vanishing entries in the working vectors. With the original algorithm, in the mean of all entries are nonzero after the second iteration, while with the improved algorithm only of the original  dimensions remain in the mean. This trend continues in subsequent iterations such that the final working dimensionality is reached more quickly with the algorithm proposed in this paper. Although using Algorithm 2 to perform the simplex projection is more expensive than just setting negative entries to zero in the orthant projection, the overhead quickly amortizes because of the boost in dimensionality reduction.

For determination of the relative speedup incorporated with both the simplex projection and the access to unit-stride arrays due to the permutation-invariance, both algorithms were implemented as C++ programs using an optimized implementation of the BLAS library for carrying out the vector operations. The employed processor was an Intel Core i7-990X. For a range of different dimensionalities, a set of vectors with varying initial sparseness were sampled. The number of the vectors for every pair of dimensionality and initial sparseness was chosen such that the processing time of the algorithms was several orders of magnitudes greater than the latency time of the operation system. Then the absolute time needed for the algorithms to compute the projections with a target sparseness of  were measured, and their ratio was taken to compute the relative speedup. The results of this experiment are depicted in Figure 3. It is evident that the maximum speedup is achieved for vectors with a dimensionality between and , and an initial sparseness greater than . For low initial sparseness, as is achieved by randomly sampled vectors, a speedup of about can be achieved for a broad spectrum of dimensionality between and .

The improvements to the original algorithm are thus not only theoretical, but also noticeable in practice. The speedup is especially useful when the projection is used as a neuronal transfer function in a classifier as proposed in Section 3, because then the computational complexity of the prediction of class membership of unknown samples can be reduced.

2.4 Function Definition and Differentiability

It is clear from Theorem 2.1 that the projection onto is unique almost everywhere. Therefore the set is a null set. However, as for example the projection is not unique for vectors where all entries are identical. In other words, for for some follows and . If a possible solution is given by with and given as stated in Remark C.1.2, as in this case and are positive. Additionally, another solution is given by which is unequal to the other solution because of . A similar argument can be used to show non-uniqueness for all . As is merely a small set, non-uniqueness is not an issue in practical applications.

The sparseness-enforcing projection operator that is restricted to non-negative solutions can thus be cast almost everywhere as a function

Exploiting reflection-invariance implies that the unrestricted variant of the projection

is well-defined, where is given as described in Lemma B. Note that computation of is a crucial prerequisite to computation of the unrestricted variant . It will be used exclusively in Section 3 because non-negativity is not necessary in the application proposed there.

If or is employed in an objective function that is to be optimized, the information whether these functions are differentiable is crucial for selecting an optimization strategy. As an example, consider once more projections onto where is a constant. It was already mentioned in Section 2 that the projection onto consists simply of zeroing out the elements that are smallest in absolute value. Let be a point and let be a permutation such that . Clearly, if then where for and for . Moreover, when then there exists a neighborhood of such that for all . With this closed-form expression, is differentiable in with gradient

, that is the identity matrix where the entries on the diagonal belonging to small absolute values of

have been zeroed out. If the requirement on is not fulfilled, then a small distortion of is sufficient to find a point in which the projection onto is differentiable.

In contrast to the projection, differentiability of and is non-trivial. A full-length discussion is given in Appendix D, and concludes that both and are differentiable almost everywhere. It is more efficient when only the product of the gradient with an arbitrary vector needs to be computed, see Corollary D

. Such an expression emerges in a natural way by application of the chain rule to an objective function where the sparseness-enforcing projection is used. In practice this weaker form is thus mostly no restriction and preferable for efficiency reasons over the more general complete gradient as given in Theorem 


The derivative of is obtained by exploiting the structure of Algorithm 1. Because the projection onto is essentially a composition of projections onto , , and , the overall gradient can be computed using the chain rule. The gradients of the intermediate projections are simple expressions and can be combined to yield one matrix for each iteration of alternating projections. Since these iteration gradients are basically sums of dyadic products, their product with an arbitrary vector can be computed by primitive vector operations. With matrix product associativity, this process can be repeated to efficiently compute the product of the gradient of with an arbitrary vector. For this, it is sufficient to record some intermediate quantities during execution of Algorithm 3, which does not add any major overhead to the algorithm itself. The gradient of the unrestricted variant can be deduced in a straightforward way from the gradient of because of their close relationship.

3 Sparse Activity and Sparse Connectivity in Supervised Learning

The sparseness-enforcing projection operator can be cast almost everywhere as vector-valued function , which is differentiable almost everywhere, see Section 2.4. This section proposes a hybrid of an auto-encoder network and a two-layer neural network, where the sparseness projection is employed as a neuronal transfer function. The proposed model is called supervised online auto-encoder (SOAE) and is intended for classification by means of a neural network that features sparse activity and sparse connectivity. Because of the analytical properties of the sparseness-enforcing projection operator, the model can be optimized end-to-end using gradient-based methods.

3.1 Architecture

Figure 4 depicts the data flow in the proposed model. There is one module for reconstruction capabilities and one module for classification capabilities. The reconstruction module, depicted on the left of Figure 4, operates by converting an input sample into an internal representation , and then computing an approximation to the original input sample. In doing so, the product of the input sample with a matrix of bases is computed, and a transfer function is applied. For sparse activity, can be chosen to be the sparseness-enforcing projection operator or the projection with respect to the pseudo-norm. This guarantees that the internal representation is sparsely populated and close to . The reconstruction is achieved like in a linear generative model, by multiplication of the matrix of bases with the internal representation. Hence the same matrix

is used for both encoding and decoding, rendering the reconstruction module symmetric, or in other words with tied weights. This approach is similar to principal component analysis

(Hotelling, 1933)

, restricted Boltzmann machines for deep auto-encoder networks

(Hinton et al., 2006) and to sparse encoding symmetric machine (Ranzato et al., 2008).

By enforcing to be sparsely populated, the sparse connectivity property holds as well. More formally, the aim is that holds for all , where is the target degree of connectivity sparseness and is the -th column of . This condition was adopted from non-negative matrix factorization with sparseness constraints (Hoyer, 2004). In the context of neural networks, the synaptic weights of individual neurons are stored in the columns of the weight matrix . The interpretation of this formal sparseness constraint is then that each neuron is only allowed to be sparsely connected with the input layer.

Figure 4: Architecture and data flow of supervised online auto-encoder (SOAE). The circle on the left, mapping from an input sample to its approximation comprises the reconstruction module. The classification module consists of the mapping from to the classification decision . The matrix of bases shall be sparsely populated to account for the sparse connectivity property. If the transfer function is set to the sparseness projection, the internal representation will be sparsely populated, fulfilling the sparse activity property.

The classification module is shown on the right-hand side of Figure 4. It computes a classification decision by feeding through a one-layer neural network. The network output is yielded through computation of the product with a matrix of weights , addition of a threshold vector and application of a transfer function . This module shares the inference of the internal representation with the reconstruction module, which can also be considered a one-layer neural network. Therefore the entire processing path from to forms a two-layer neural network (Rumelhart et al., 1986), where stores the synaptic weights of the hidden layer, and and are the parameters of the output layer.

The input sample shall be approximated by , and the target vector for classification shall be approximated by . This is achieved by optimization of the parameters of SOAE, that is the quantities , and . The goodness of the approximation is estimated using a differentiable similarity measure , and the approximation is assessed by another similarity measure . For minimizing the deviation in both approximations, the objective function

shall be optimized, where controls the trade-off between reconstruction and classification capabilities. To incorporate sparse connectivity, feasible solutions are restricted to fulfill for all . If , then SOAE is identical to a symmetric auto-encoder network with sparse activity and sparse connectivity. In the case of , SOAE forms a two-layer neural network for classification with a sparsely connected hidden layer and where the activity in the hidden layer is sparse. The parameter can also be used to blend continuously between these two extremes. Note that only depends on but not on or , but depends on , and . Hence and are only relevant when , whereas is essential for all choices of .

An appropriate choice for is the correlation coefficient (see for example Rodgers and Nicewander, 1988), because it is normed to values in the interval

, invariant to affine-linear transformations, and differentiable. If

is set to , then a model that is invariant to the concrete scaling and shifting of the occurring quantities can be yielded. This follows because is also invariant to such transformations, see Corollary C.1.2. The similarity measure for classification capabilities is chosen to be the cross-entropy error function (Bishop, 1995), which was shown empirically by Simard et al. (2003) to induce better classification capabilities than the mean squared error function. The softmax transfer function (Bishop, 1995) is used as transfer function of the output layer. It provides a natural pairing together with the cross-entropy error function (Dunne and Campbell, 1997) and supports multi-class classification.

3.2 Learning Algorithm

The proposed optimization algorithm for minimization of the objective function is projected gradient descent (Bertsekas, 1999)

. Here, each update to the degrees of freedom is followed by application of the sparseness projection to the columns of

to enforce sparse connectivity. There are theoretical results on the convergence of projected gradient methods when projections are carried out onto convex sets (Bertsekas, 1999)

, but here the target set for projection is non-convex. Nevertheless, the experiments described below show that projected gradient descent is an adequate heuristic in the situation of the SOAE framework to tune the network parameters. For completeness, the gradients of

with respect to the network parameters are given in Appendix E. Update steps are carried out after every presentation of a pair of an input sample and associated target vector. This online learning procedure results in faster learning and improves generalization capabilities over batch learning (Wilson and Martinez, 2003; Bottou and LeCun, 2004).

A learning set with samples from and associated target vectors from as one-of--codes is input to the algorithm. The dimensionality of the internal representation and the target degree of sparseness with respect to the connectivity are parameters of the algorithm. Sparseness of connectivity increases for larger , as Hoyer’s sparseness measure is employed in the definition of the set of feasible solutions.

Two possible choices for the hidden layer’s transfer function to achieve sparse activity were discussed in this paper. One possibility is to carry out the projection with respect to the pseudo-norm. The more sophisticated method is to use the unrestricted sparseness-enforcing projection operator with respect to Hoyer’s sparseness measure , which can be carried out by Algorithm 3. In both cases, a target degree for sparse activity is a parameter of the learning algorithm. In case of the projection, this sparseness degree is denoted by , and sparseness increases with smaller values of it. For the projection, is used, where larger values indicate more sparse activity.

Initialization of the columns of

is achieved by selecting a random subset of the learning set, similar to the initialization of radial basis function networks

(Bishop, 1995). This ensures significant activity of the hidden layer from the very start, resulting in strong gradients and therefore reducing training time. The parameters of the output layer, that is and

, are initialized by sampling from a zero-mean Gaussian distribution with a standard deviation of


In every epoch, a randomly selected subset of samples and associated target vectors from the learning set is used for stochastic gradient descent to update

, and . The results from Appendix E can be used to efficiently compute the gradient of the objective function. There, the gradient for the transfer function only emerges as a product with a vector. The gradient for the projection is trivial and was given as an example in Section 2.4. If is Hoyer’s sparseness-enforcing projection operator, it is possible to exploit that only the product of the gradient with a vector is needed. In this case, it is more efficient to compute the result of the multiplication implicitly using Corollary D and thus avoid the computation of the entire gradient of .

After every epoch, a sparseness projection is applied to the columns of . This guarantees that holds for all , and therefore the sparse connectivity property is fulfilled. The trade-off variable which controls the weight of the reconstruction and the classification term is adjusted according to , where denotes the number of the current epoch. Thus starts at zero, increases slowly and asymptotically reaches one. The emphasis at the beginning of the optimization is thus on reconstruction capabilities. Subsequently, classification capabilities are incorporated slowly, and in the final phase of training classification capabilities exclusively are optimized. This continuous variant of unsupervised pre-training (Hinton et al., 2006) leads to parameters in the vicinity of a good minimizer for classification capabilities before classification is preferred over reconstruction through the trade-off parameter . Compared to the choice this strategy helps to stabilize the trajectory in parameter space and makes the objective function values settle down more quickly, such that the termination criterion is satisfied earlier.

3.3 Description of Experiments

To assess the classification capabilities and the impact of sparse activity and sparse connectivity, the MNIST database of handwritten digits (LeCun and Cortes, 1998) was employed. It is a popular benchmark data set for classification algorithms, and numerous results with respect to this data set are reported in the literature. The database consists of samples, divided into a learning set of samples and an evaluation set of samples. Each sample represents a digit of size pixels and has a class label from associated with it. Therefore the input and output dimensionalities are and , respectively. The classification error is given in percent of all evaluation samples, hence corresponds to a single misclassified digit.

For generation of the original data set, the placement of the digits has been achieved based on their barycenter (LeCun and Cortes, 1998). Because of sampling and rounding errors, the localization uncertainty can hence be assumed to be less than one pixel in both directions. To account for this uncertainty, the learning set was augmented by jittering each sample in each of eight possible directions by one pixel, yielding samples for learning in total. The evaluation set was left unchanged to yield results that can be compared to the literature. As noted by Hinton et al. (2006), the learning problem is no more permutation-invariant due to the jittering, as information on the neighborhood of the pixels is implicitly incorporated in the learning set.

However, classification results improve dramatically when such prior knowledge is used. This was demonstrated by Schölkopf (1997)

using the virtual support vector method, which improved a support vector machine with polynomial kernel of degree five from an error of

to by jittering the support vectors by one pixel in four principal directions. This result was extended by DeCoste and Schölkopf (2002), where a support vector machine with a polynomial kernel of degree nine was improved from an error of to by jittering in all possible eight directions. Further improvements can be achieved by generating artificial training samples using elastic distortions (Simard et al., 2003). This reduced the error of a two-layer neural network with  hidden units to , compared to the error yielded when training on samples created by affine distortions. Very big and very deep neural networks possess a large number of adaptable weights. In conjunction with elastic and affine distortions such neural networks can yield errors as low as (Cireşan et al., 2010). The current record error of

is held by an approach that combines distorted samples with a committee of convolutional neural networks

(Cireşan et al., 2012). This is an architecture that has been optimized exclusively for input data that represents images, that is where the neighborhood of the pixels is hard-wired in the classifier. To allow for a plain evaluation that does not depend on additional parameters for creating artificial samples, the jittered learning set with samples is used throughout this paper.

The experimental methodology was as follows. The number of hidden units was chosen to be in all experiments that are described below. This is an increased number compared to the  hidden units employed by Simard et al. (2003), but promises to yield better results when an adequate number of learning samples is used. As all tested learning algorithms are essentially gradient descent methods, an initial step size had to be chosen. For each candidate step size, five runs of a two-fold cross validation were carried out on the learning set. Then, for each step size the median of the ten resulting classification errors was computed. The winning step size was then determined to be the one that achieved a minimum median of classification errors.

In every epoch, samples were randomly chosen from the learning set and presented to the network. This number of samples was chosen as it is -th of the jittered learning set. The step size was multiplicatively annealed using a factor of after every epoch. Optimization was terminated once the relative change in the objective function became very small and no more significant progress on the learning set could be observed. The resulting classifiers were then applied to the evaluation set, and misclassifications were counted.

3.4 Experimental Results

Two variants of the supervised online auto-encoder architecture as proposed in this section were trained on the augmented learning set. In both variants, the target degree of sparse connectivity was set to . This choice was made because  of all samples in the learning set possess a sparseness which is less than . Therefore, the resulting bases are forced to be truly sparsely connected compared to the sparseness of the digits.

The first variant is denoted by SOAE-. Here, the sparseness-enforcing projection operator was used as transfer function in the hidden layer. Target degrees of sparse activity with respect to Hoyer’s sparseness measure were chosen from the interval in steps of size . This variant was then trained on the jittered learning set using the method described in Section 3.2. For every value of , the resulting sparseness of activity was measured after training using the pseudo-norm. For this, each sample of the learning set was presented to the networks, and the number of active units in the hidden layer was counted. Figure 5 shows the resulting mean value and standard deviation of sparse activity. If is chosen, then in the mean about  of the total  hidden units are active upon presentation of a sample from the learning set. For only one hundred units are active at any one time, and for there are only eleven active units. The standard deviation of the activity decreases when sparseness increases, hence the mapping from to the resulting number of active units becomes more accurate.

Figure 5: Resulting amount of nonzero entries in an internal representation with  entries, depending on the target degree of sparseness for activity with respect to Hoyer’s sparseness measure . For low values of , about of the entries are nonzero, whereas for very high sparseness degrees only of the entries do not vanish. The error bars indicate  one standard deviation distance from the mean value. Standard deviation shrinks with increasing sparseness degree, making the mapping more accurate.

The second variant, denoted SOAE-, differs from SOAE- in that the projection with respect to the pseudo-norm as transfer function was used. The target sparseness of activity is given by a parameter , which controls the exact number of units that are allowed to be active at any one time. For the experiments, the values for were chosen to match the mean activities from the SOAE- experiments. This way the results of both variants can be compared based on a unified value of activity sparseness. The results are depicted in Figure 6. Usage of the projection consequently outperforms the projection for all sparseness degrees. Even for high sparseness of activity, that is when only about ten percent of the units are allowed to be active at any one time, good classification capabilities can be obtained with SOAE-. For , the classification results of SOAE- reach an optimum. SOAE- is more robust, as classification capabilities first begin to collapse when sparseness is below , whereas SOAE- starts to degenerate when sparseness falls below . For , roughly translating to between and activity, about equal classification performance is achieved using SOAE-.

Figure 6: Resulting classification error on the MNIST evaluation set for the supervised online auto-encoder network, in dependence of sparseness of activity in the hidden layer. The projection onto an pseudo-norm constraint for variant SOAE- and the projection onto a constraint determined by sparseness measure for variant SOAE- were used as transfer functions. The error bars indicate  one standard deviation difference from the mean.

It can thus be concluded that using the sparseness-enforcing projection operator as described in this paper yields better results than when the simple projection is used to achieve sparse activity. To assess the benefit more precisely and to investigate the effect of individual factors, several comparative experiments have been carried out. A summary of these experiments and their outcome is given in Table 1. The variants SOAE- and SOAE- denote the entirety of the respective experiments where sparseness of activity lies in the intervals described above, that is and , respectively. Using these intervals, SOAE- and SOAE- achieved a median error of and on the evaluation set, respectively. Variant SOAE--conn is essentially equal to SOAE-, except for sparse connectivity not being incorporated. Sparseness of activity here was also chosen to be , which resulted in about equal classification results over the entire range. Dropping of sparse connectivity increases misclassifications, for the median error of SOAE--conn is and thereby greater than the median error of SOAE-.

Approach Sparse Sparse Result of Evaluation
Connectivity Activity Shapiro-Wilk Test Error [%]
SOAE--conn none
MLP-OBD none
MLP-random none none
MLP-samples none none
MLP-SCFC none none
Table 1: Overview of comparative experiments. The second and third columns indicate whether sparse connectivity or sparse activity was incorporated, respectively. The fourth column reports the result of a statistical test for normality, which is interpreted in Section 3.5. The final column gives the median one standard deviation of the achieved classification error on the MNIST evaluation set. The results for each experiment were trimmed to gain a sample of size , allowing for statistical robust estimates.

The other five approaches included in the comparison are multi-layer perceptrons (MLPs) with the same topology and dynamics as the classification module of supervised online auto-encoder, with two exceptions. First, the transfer function of the hidden layer

was set to a hyperbolic tangent, thus not including explicit sparse activity. Second, in all but one experiment sparse connectivity was either not incorporated, or achieved through other means than by performing a projection after each learning epoch. Besides the variation in sparseness of connectivity, the experiments differ in the initialization of the network parameters.

For each variant, runs were carried out and the resulting classifiers were applied to the evaluation set to compute the classification error. Then, the best four and the worst four results were discarded and not included in further analysis. Hence a random sample of size was achieved, where of the original data were trimmed away. This procedure was also applied to the results of SOAE-, SOAE--conn, and SOAE-, to obtain a total of eight random samples of equal size for comparison with another.

The most basic variant, denoted the baseline in this discussion, is MLP-random, where all network parameters were initialized randomly. This achieved a median error of on the evaluation set, being considerably worse than SOAE-. For variant MLP-samples, the hidden layer was initialized by replication of randomly chosen samples from the learning set. This did decrease the overall learning time. However, the median classification error was slightly worse with compared to MLP-random.

For variant MLP-SCFC, the network parameters were initialized in an unsupervised manner using the sparse coding for fast classification (SCFC) algorithm (Thom et al., 2011a). This method is a precursor to the SOAE proposed in this paper. It also features sparse connectivity and sparse activity but differs in some essential parts. First, sparseness of activity is achieved through a latent variable that stores the optimal sparse code words of all samples simultaneously. Using this matrix of code words, the activity of individual units was enforced to be sparse over time on the entire learning set. SOAE achieves sparseness over space, as for each sample only a pre-defined fraction of units is allowed to be active at any one time. A second difference is that sparse activity is achieved only indirectly by approximation of the latent matrix of code words with a feed-forward representation. With SOAE, sparseness of activity is guaranteed by construction. MLP-SCFC achieved a median classification error of on the MNIST evaluation set, rendering it slightly worse than MLP-random and equivalent to MLP-samples.

The first experiment that incorporates only sparse connectivity is SMLP-SCFC. Initialization was done as for MLP-SCFC, but during training sparseness of connectivity was yielded by application of the sparseness-enforcing projection operator to the weights of the hidden layer after every learning epoch. Hence the sparseness gained from unsupervised initialization was retained. MLP-SCFC features sparse connectivity only after initialization, but loses this property when training proceeds. With this slight modification, the median error of SMLP-SCFC decreases to , which is significantly better than the baseline result.

The effect of better generalization due to sparse connectivity has also been observed by LeCun et al. (1990)

in the context of convolutional neural networks. It can be explained by the bias-variance decomposition of the generalization error

(Geman et al., 1992). When the effective number of the degrees of freedom is constrained, overfitting will be less likely and hence classifiers produce better results on average. The same argument can be applied to SOAE-, where additional sparse activity further improves classification results.

The last variant is called MLP-OBD. Here, the optimal brain damage (OBD) algorithm (LeCun et al., 1990) was used to prune synaptic connections in the hidden layer that are irrelevant for the computation of the classification decision of the network. The parameters of the network were first initialized randomly and then optimized on the learning set. Then the impact for each synaptic connection on the objective function was estimated using the Taylor series of the objective function, where a diagonal approximation of the Hessian was employed and terms of cubic or higher order were neglected. Using this information, the number of connections was halved by setting the weight of connections with low impact to zero. The network was then retrained with weights of removed connections kept at zero. This procedure was repeated until a target percentage of active synaptic connections in the hidden layer was achieved. For the results reported here, was chosen as this reflects the sparse connectivity of the other approaches best. MLP-OBD achieved a median classification error of , which is comparable to the baseline result.

3.5 Statistical Analysis and Conclusions

A statistical analysis was carried out to assess the significance of the differences in the performance of the eight algorithms. The procedure follows the proposals of Pizarro et al. (2002) and Demšar (2006) for hypothesis testing, and is concluded by effect size estimation as proposed by Grissom (1994) and Acion et al. (2006). For each algorithm, a sample of size was available, allowing for robust analysis results.

First, all results were tested for normality using the test developed by Shapiro and Wilk (1965)

. The resulting test statistics

and -values are given in Table 1. As all

-values are large, it cannot be rejected that the samples came from normally distributed populations. Thus normality is assumed in the remainder of this discussion. Next, the test proposed by

Levene (1960) was applied to determine whether equality of variances of the groups holds. This resulted in a test statistic with and degrees of freedom, and therefore a -value of . Hence the hypothesis that all group variances are equal can be rejected with very high significance. Consequently, parametric omnibus and post-hoc tests cannot be applied, as they require the groups to have equal variance.

As an alternative, the nonparametric test by Kruskal and Wallis (1952) which is based on rank information was employed to test whether all algorithms produced classifiers with equal classification errors in the mean. The test statistic was with degrees of freedom, and the -value was less than

. There is hence a statistically significant difference in the mean classification results. To locate this deviation, a critical difference for comparing the mean ranks of the algorithms was computed. A Tukey-Kramer type modification applied to Dunn’s procedure yields this critical difference, which is less conservative than Nemenyi’s procedure for the Kruskal-Wallis test

(Hochberg and Tamhane, 1987). Note that this approach is nevertheless similar to the post-hoc procedure proposed by Demšar (2006) for paired observations, such that the diagrams proposed there can be adapted to the case for unpaired observations. The result is depicted in Figure 7, where the critical difference for statistical significance at the level is given. This test induces a highly significant partitioning of the eight algorithms, namely three groups , and given by

This partition in turn induces an equivalence relation. Statistical equivalence is hence unambiguous and well-defined at . Moreover, the -value for this partition is . If the significance level would have been set lower than this, then groups and would blend together.

Figure 7: Diagram for multiple comparison of algorithms following Demšar (2006). For each algorithm, the mean rank was computed during the Kruskal-Wallis test. Then, a critical difference (CD) was computed at the significance level. Two algorithms produce classification results that are statistically not equal if the difference between their mean ranks is greater than the critical difference. This induced three groups of algorithms that produced statistically equivalent results, which are marked with black bars.

To assess the benefit when an algorithm from one group is chosen over an algorithm from another group, the probability of superior experiment outcome was estimated

(Grissom, 1994; Acion et al., 2006). For this, the classification errors were pooled with respect to membership in the three groups. It was then tested whether these pooled results still come from normal distributions. As group is a singleton, this is trivially fulfilled with the result from Table 1. For group , the Shapiro-Wilk test statistic was and the -value was . Group achieved a test statistic of and a -value of . If a standard significance level of is chosen, then and can be assumed to be normally distributed also.


be the random variable modeling the classification results of the algorithms from group

. It is assumed that is normally distributed with unknown mean and unknown variance for all . Then is clearly normally distributed also for two groups . Therefore, the probability

that one algorithm produces a better classifier than another could be computed from the Gaussian error function if the group means and variances were known. However, using Rao-Blackwell theory a minimum variance unbiased estimator

of this probability can be computed easily (Downton, 1973). Evaluation of the expression for shows that can be estimated by , can be estimated by , and can be estimated by . Therefore, the effect of choosing SOAE- over any of the seven other algorithms is dramatic (Grissom, 1994).

These results can be interpreted as follows. When neither sparse activity nor sparse connectivity is incorporated, then the worst classification results are obtained regardless of the initialization of the network parameters. The exception is MLP-OBD which incorporates sparse connectivity, although, as its name says, in a destructive way. Once a synaptic connection has been removed, it cannot be recovered, as the measure for relevance of LeCun et al. (1990) vanishes for synaptic connections of zero strength. The statistics for SMLP-SCFC shows that when sparse connectivity is obtained using the sparseness-enforcing projection operator, then superior results can be achieved. Because of the nature of projected gradient descent, it is possible here to restore deleted connections if it helps to decrease the classification error during learning. For SOAE--conn only sparse activity was used, and classification results were statistically equivalent to SMLP-SCFC.

Therefore, using either sparse activity or sparse connectivity improves classification capabilities. When both are used, then results improve even more as variant SOAE- shows. This does not hold for SOAE- however, where the projection was used as transfer function. As Hoyer’s sparseness measure and the according projection possess desirable analytical properties, they can be considered smooth approximations to the pseudo-norm. It is this smoothness which seems to produce this benefit in practice.

4 Related Work

This section reviews work related with the contents of this paper. First, the theoretical foundations of the sparseness-enforcing projection operator are discussed. Next, its application as neuronal transfer function to achieve sparse activity in a classification scenario is put in context with alternative approaches, and possible advantages of sparse connectivity are described.

4.1 Sparseness-Enforcing Projection Operator

The first major part of this paper dealt with improvements to the work of Hoyer (2004) and Theis et al. (2005). Here, an algorithm for the sparseness-enforcing projection with respect to Hoyer’s sparseness measure was proposed. The technical proof of correctness is given in Appendix C. The set that should be projected onto is an intersection of a simplex and a hypercircle , which is a hypersphere lying in a hyperplane. The overall procedure can be described as performing alternating projections onto and certain subsets of . This approach is common for handling projections onto intersections of individual sets. For example, von Neumann (1950) proposed essentially the same idea when the investigated sets are closed subspaces, and has shown that this converges to a solution. A similar approach can be carried out for intersections of closed, convex cones (Dykstra, 1983), which can be generalized to translated cones that can be used to approximate any convex set (Dykstra and Boyle, 1987). For these alternating methods, it is only necessary to know how projections onto individual members of the intersection can be achieved.

Although these methods exhibit great generality, they have two severe drawbacks in the scenario of this paper. First, the target set for projection must be an intersection of convex sets. The scaled canonical simplex is clearly convex, but the hypercircle is non-convex if it contains more than one point. The condition that generates cannot easily be weakened to achieve convexity. If the original hypersphere were replaced with a closed ball, then would be convex. But this changes the meaning of the problem dramatically, as now virtually any sparseness below the original target degree of sparseness can be obtained. This is because when the target norm is fixed, the sparseness measure decreases whenever the target norm decreases. In geometric terms, the method proposed in this paper performs a projection from within a circle onto its boundary to increase the sparseness of the working vector. This argument is given in more detail in Figure 11 and the proof of Lemma C.36.

The second drawback of the general methods for projecting onto intersections is that a solution is only achieved asymptotically, even when the convexity requirements are fulfilled. Due to the special structure of and , the number of alternating projections that have to be carried out to find a solution using Algorithm 3 is bounded from above by the problem dimensionality. Thus an exact projection is always found in finite time. Furthermore, the solution is guaranteed to be found in time that is at most quadratic in problem dimensionality.

A crucial point is the computation of the projection onto and certain subsets of . Due to the nature of the norm, the latter is straightforward. For the former, efficient algorithms have been proposed recently (Duchi et al., 2008; Chen and Ye, 2011). When only independent solutions are required, the projection of a point onto a scaled canonical simplex of norm can also be carried out in linear time (Liu and Ye, 2009), without having to sort the vector that is to be projected. This can be achieved by showing that the separator for performing the simplex projection is the unique zero of the monotonically decreasing function . The zero of this function can be found efficiently using the bisection method, and exploiting the special structure of the occurring expressions (Liu and Ye, 2009).

In the context of this paper an explicit closed-form expression for is preferable as it permits additional insight into the properties of the projected point. The major part in proving the correctness of Algorithm 1 is the interconnection between and , that is that the final solution has zero entries at the according positions in the working vector and thus a chain monotonically decreasing in pseudo-norm is achieved. This result is established through Lemma C.2.2, which characterizes projections onto certain faces of a simplex, Corollary 10 and their application in Lemma C.3.

Analysis of the theoretical properties of the sparseness-enforcing projection is concluded with its differentiability in Appendix D. The idea is to exploit the finiteness of the projection sequence and to apply the chain rule of differential calculus. It is necessary to show that the projection chain is robust in a neighborhood of the argument. This reduces analysis to individual projection steps which have already been studied in the literature. For example, the projection onto a closed, convex set is guaranteed to be differentiable almost everywhere (Hiriart-Urruty, 1982). Here non-convexity of is not an issue, as the only critical point is its barycenter. For the simplex , a characterization of critical points is given with Lemma D and Lemma D, and it is shown that the expression for the projection onto is invariant to local changes. An explicit expression for construction of the gradient of the sparseness-enforcing projection operator is given in Theorem D. In Corollary D it is shown that the computation of the product of the gradient with an arbitrary vector can be achieved efficiently by exploiting sparseness and the special structure of the gradient.

Similar approaches for sparseness projections are discussed in the following. The iterative hard thresholding algorithm is a gradient descent algorithm, where a projection onto an pseudo-norm constraint is performed (Blumensath and Davies, 2009). Its application lies in compressed sensing, where a linear generative model is used to infer a sparse representation for a given observation. Sparseness here acts as regularizer which is necessary because observations are sampled below the Nyquist rate. In spite of the simplicity of the method, it can be shown that it achieves a good approximation to the optimal solution of this NP-hard problem (Blumensath and Davies, 2009).

Closely related with the work of this paper is the generalization of Hoyer’s sparseness measure by Theis and Tanaka (2006). Here, the norm constraint is replaced with a generalized pseudo-norm constraint, such that the sparseness measure becomes . For , Hoyer’s sparseness measure up to a constant normalization is obtained. When converges decreasingly to zero, then converges point-wise to the pseudo-norm. Hence for small values of a more natural sparseness measure is obtained. Theis and Tanaka (2006) also proposed an extension of Hoyer’s projection algorithm. It is essentially von Neumann’s alternating projection method, where closed subspaces have been replaced by "spheres" that are induced by pseudo-norms. Note that these sets are non-convex when , such that convergence is not guaranteed. Further, no closed-form solution for the projection onto an "-sphere" is known for , such that numerical methods have to be employed.

A problem where similar projections are employed is to minimize a convex function subject to group sparseness (see for example Friedman et al., 2010). In this context, mixed norm balls are of particular interest (Sra, 2012). For a matrix , the mixed norm is defined as the norm of the norms of the columns of , that is . Here, can be interpreted to be a data point with entries partitioned into groups. When , then the projection onto a simplex can be generalized directly for (van den Berg et al., 2008) and for (Quattoni et al., 2009). The case when and is more difficult, but can be solved as well (Liu and Ye, 2010; Sra, 2012).

The last problem discussed here is the elastic net criterion (Zou and Hastie, 2005), which is a constraint on the sum of an norm and an norm. The feasible set can be written as the convex set , where control the shape of . Note that in only the sum of two norms is considered, whereas the non-convex set consists of the intersection of two different constraints. Therefore, the elastic net induces a different notion of sparseness than Hoyer’s sparseness measure does. As is the case for mixed norm balls, the projection onto a simplex can be generalized to achieve projections onto (Mairal et al., 2010).

4.2 Supervised Online Auto-Encoder

The sparseness-enforcing projection operator with respect to Hoyer’s sparseness measure and the projection onto an pseudo-norm constraint are differentiable almost everywhere. Thus they are suitable for gradient-based optimization algorithms. In Section 3, they were used as transfer functions in a hybrid of an auto-encoder network and a two-layer neural network to infer a sparse internal representation. This representation was subsequently employed to approximate the input sample and to compute a classification decision. In addition, the matrix of bases which was used to compute the internal representation was enforced to be sparsely populated by application of the sparseness projection after each learning epoch. Hence the supervised online auto-encoder proposed in this paper features both sparse activity and sparse connectivity.

These two key properties have also been investigated and exploited in the context of auto-associative memories for binary inputs. If the entries of the training patterns are sparsely populated, the weight matrix of the memory will be sparsely populated as well after training if Hebbian-like learning rules are used

(Kohonen, 1972). The assumption of sparsely coded inputs also results in increased completion capacity and noise resistance of the associative memory (Palm, 1980). If the input data is not sparse inherently, feature detectors can perform a sparsification prior to the actual processing through the memory (Baum et al., 1988).

A purely generative model that also possesses these two key properties is non-negative matrix factorization with sparseness constraints (Hoyer, 2004). This is an extension to plain non-negative matrix factorization (Paatero and Tapper, 1994) which was shown to achieve sparse connectivity on certain data sets (Lee and Seung, 1999). However, there are data sets on which this does not work (Li et al., 2001; Hoyer, 2004). Although Hoyer’s model makes sparseness easily controllable by explicit constraints, it is not inherently suited to classification tasks. An extension intended to incorporate class membership information to increase discriminative capabilities was proposed by Heiler and Schnörr (2006). In their approach, an additional constraint was added ensuring that every internal representation is close to the mean of all internal representations that belong to the same class. In other words, the method can be interpreted as supervised clustering, with the number of clusters equal to the number of classes. However, there is no guarantee that a distribution of internal representations exists such that both the reproduction error is minimized and the internal representations can be arranged in such a pattern. Unfortunately, Heiler and Schnörr (2006) used only a subset of a small data set for handwritten digit recognition to evaluate their approach.

A precursor to the supervised online auto-encoder was proposed by Thom et al. (2011a)

. There, inference of sparse internal representations was achieved by fitting a one-layer neural network to approximate a latent variable of optimal sparse representations. The transfer function used for this approximation was a hyperbolic tangent raised to an odd power greater or equal to three. This resulted in a depression of activities with small magnitude, favoring sparseness of the result. Similar techniques to achieve a shrinkage-like effect for increasing sparseness of activity in a neural network were used by

Gregor and LeCun (2010) and Glorot et al. (2011). Information processing is here purely local, that is a scalar function is evaluated entrywise on a vector, and thus no information is interchanged among individual entries.

The use of non-local shrinkage to reduce Gaussian noise in sparse coding has already been described by Hyvärinen et al. (1999). Here, a maximum likelihood estimate with only weak assumptions yields a shrinkage operation, which can be conceived as projection onto a scaled canonical simplex. In the use case of object recognition, a hard shrinkage was also employed to de-noise filter responses (Mutch and Lowe, 2006). Whenever a best approximation from a permutation-invariant set is used, a shrinkage-like operation must be employed. Using a projection operator as neural transfer function is hence a natural extension of these ideas. When the projection is sufficiently smooth, the entire model can be tuned end-to-end using gradient methods to achieve an auto-encoder or a classifier.

The second building block from Thom et al. (2011a) that was incorporated into supervised online auto-encoder is the architectural concept for classification. It is well-known that two layers in a neural network are sufficient to approximate any continuous function on a compactum with arbitrary precision (Cybenko, 1989; Funahashi, 1989; Hornik et al., 1989). Similar architectures have also been proposed for classification in combination with sparse coding of the inputs. However, sparse connectivity was not considered in this context. Bradley and Bagnell (2009)

used the Kullback-Leibler divergence as implicit sparseness penalty term and combined this with the backpropagation algorithm to yield a classifier that achieved a

error rate on the MNIST evaluation set. The Kullback-Leibler divergence was chosen to replace the usual norm penalty term, as it is smoother than the latter and therefore sparsely coded internal representations are more stable subject to subtle changes of the input. A related technique is supervised dictionary learning by Mairal et al. (2009), where the objective function is an additive combination of a classification error term, a term for the reproduction error, and an norm constraint. Inference of sparse internal representations is achieved through solving an optimization problem. Such procedures are time-consuming and greatly increase the computational complexity of classification. With this approach, a classification error of on the MNIST evaluation set was achieved. These two approaches used the original MNIST learning set without jittering the digits and can thus be considered permutation-invariant. Augmentation of the learning set with virtual samples would have contributed to improve classification performance, as demonstrated by Schölkopf (1997).

Finally consider once more the sparse connectivity property, which is mostly neglected in the literature in favor of sparse activity. It was shown in this paper that sparse connectivity helps to improve generalization capabilities. In practice, this property can also be used to reduce the computational complexity of classification by one order of magnitude (Thom et al., 2011b). This results from exploiting sparseness and using sparse matrix-vector multiplication algorithms to infer the internal representation, which is the major computational burden in class membership prediction. It was shown in this paper and by Thom et al. (2011b) that a small number of nonzero entries in the weight matrix of the hidden layer is sufficient for achieving good classification results. Furthermore, the additional savings in required storage capacity and bandwidth allow using platforms with modest computational power for practical implementations. Sparseness is therefore an elementary concept of efficiency in artificial processing systems.

5 Conclusions

Without sparseness in their brains, higher mammals probably would not have developed to viable life-forms. This important concept of efficiency was discovered by neuroscientists, and practical benefit was obtained by the engineers of artificial information processing systems. This paper studied Hoyer’s sparseness measure , and in particular the projection of arbitrary vectors onto sets where attains a constant value. A simple yet efficient algorithm for computing this sparseness-enforcing projection operator was proposed in this paper, and its correctness was proved. In addition, it was demonstrated that the proposed algorithm is superior in run-time to Hoyer’s original algorithm. The analysis of the theoretical properties of this projection was concluded by showing it is differentiable almost everywhere.

As projections onto constraints are well-understood, they constitute the ideal tool for building systems that can benefit from sparseness constraints. An original use case was introduced in this paper. Here, the projection was implemented as neuronal transfer function, yielding a differentiable closed-form expression for inference of sparse code words. Besides this sparse activity, the connectivity in this system was also forced to be sparse by performing the projection after the presentation of learning examples. Because of its smoothness, the entire system can be optimized end-to-end by gradient-based methods, yielding a classification architecture exhibiting true sparse information processing.

This supervised online auto-encoder was applied on a benchmark data set for pattern recognition. Because sparseness constraints reduce the amount of feasible solutions, it is not clear in the first place whether the same performance can be achieved at all. However, when the target degree of sparseness of the activity is in a reasonable range, classification results are not only equivalent but superior to classical non-sparse approaches. This result is supported by statistical evaluation showing that this performance increase is not merely coincidental, but statistically significant. Therefore, sparseness can be seen as regularizer that offers the potential to improve artificial systems in the same way it seems to improve biological systems.

The authors wish to thank Patrik O. Hoyer and Xiaojing Ye for sharing the source code of their algorithms. The authors are also grateful to the anonymous reviewers for their valuable comments and feedback. This work was supported by Daimler AG, Germany.

Appendix A Notation and Prerequisites

This appendix fixes the notation and provides prerequisites for the following appendices. denotes the natural numbers including zero, the real numbers and the non-negative real numbers. is the -dimensional Euclidean space with canonical basis , and denotes the vector where all entries are identical to unity. For all other vectors, a subscript denotes the corresponding entry of the vector, that is for . The amount of nonzero entries in a vector is given by the pseudo-norm, . and denote the Manhattan norm and Euclidean norm, respectively. denotes the canonical dot product in the Euclidean space. Given a vector , denotes the square matrix with on its main diagonal and zero entries at all other positions, and denotes the Hadamard product or entrywise product for vectors and . When and are square matrices, then denotes the block diagonal matrix with the blocks given by and . is the symmetric group, and denotes the permutation matrix for . For a set , denotes its complement in the universal set , where is clear from the context. The power set of is denoted by . If , then denotes its boundary in the topological sense. The sign function is denoted by . A list of symbols that are frequently used throughout the paper is given in Table 2.

Symbol and Definition Meaning
(see Section 1.1) Sparseness measure by Hoyer (2004)
and (see Section 2.4) Sparseness projection cast as function
Problem dimensionality
Canonical basis of
Vector where all entries are one
Target or Manhattan norm
Target or Euclidean norm
(see Section 1.1) Target set for sparseness projection
Target set for non-negative sparseness projection
Short for the non-negative target set
Target hyperplane
Target hypersphere
Target hypercircle
Scaled canonical simplex
Barycenter of and
Squared radius of
Index set of nonzero entries
Working dimensionality
Points in where certain coordinates vanish
Face of simplex
Barycenter of and