The NMF is a powerful tool used in clustering, image processing, text mining, and so on. Its importance grew in the last decade due to its efficacy into extracting meaningful and easily interpretable features out of the data. For example, in the clustering problem of points into a dimensional space, the processed data can be naturally viewed as centroids of the clusters, or in its application to text mining, the NMF output clearly points to the common topics touched by the input documents. In this paper the focus is on the applications of NMF to the analysis and decomposition of images, as shown in the article of Lee & Seung , where they processed a set of faces and the algorithm automatically recognized their principal features like eyebrows, lips, noses, etc.
A serious drawbacks of this method is that NMF can’t recognize the same objects or parts of them if they’re located in different places on multiple images, or when they’re rotated or deformed. In other words, NMF is not invariant under space transformations, so the input data must always be pre-calibrated and adjusted.
One possible solution may be to add to the dataset a lot of copy of the same image, each time stretched, rotated and shifted in different ways, in order to make the NMF recognize the parts of an image even if they’re in different positions and with different shapes, but this leads to an huge rise of input data and of redundancy in the solution.
Some authors have suggested to set some standard transformations of the images (such as translations or symmetry) and to look for the features we want to obtain, along with additional parameters that indicate for each transformation of each feature if they’re present into the original images. This rises the number of the problem variables by a factor that’s usually larger or equal to the number of pixels in a picture, like in  and , making the algorithm complexity go up by at least the same factor.
Here is presented a way to attack the problem of the translations, keeping the framework of NMF and the natural the graphical property of its output to represent the wanted parts of images, and bounding the rise in data weight and computational cost with the number of effective components we want to find and a logarithmic factor.
In the first chapter we review the original NMF problem, and we’ll discuss why it’s applicable to image processing. On the second chapter, we introduce the tools and notation needed to state the actual problem we want to solve. On the third chapter, we describe the algorithms used, and derive the asymptotic computational cost. On the fourth chapter we present some experiments on hand-made images, and on the conclusions we’ll talk about possible improvements.
1 NMF and image processing
1.1 Nonnegative Matrix Factorization
Given a data matrix and a natural number , the NMF problem requires to find the matrices that satisfy
where we used the Frobenius norm, defined as
A natural interpretation of NMF derives from the observation that, given any column of , a good solution to the problem finds an approximation of it through a combination of
nonnegative vectors, the columns of, with nonnegative coefficients stored in a row of . This means that the problem is equivalent to find a nonnegative set of vectors that approximately generate, through nonnegative coefficient, all the columns of (the minimum parameter that satisfy such conditions is often referred to as the nonnegative rank of the matrix ).
Usually, is much smaller then the other dimensions since the NMF is often used as a low-rank decomposition algorithm, and the resulting columns of , called features or components, have a meaningful representation as characteristics or parts of the original data, that are the columns of . Moreover, a large value of implies a large set of exact solutions for the exact NMF problem, and it translates into a lot of local minima into the minimization problem, that leads to inaccuracy on the algorithmic part, and ambiguity in the interpretation of solutions.
An other feature that is usually required to the input data is the sparsity, since it is proved that can improve the quality and understandability of the solution, along with gaining uniqueness properties (for further studies, see , that proposes a preprocessing to improve the sparseness of ).
A common way to take advantage form the non-uniqueness of the solution is to normalize rows and columns of , through a positive diagonal matrix of dimensions . In fact, given any pair , then , so the matrices are still nonnegative, and this transformation doesn’t change the error we want to minimize. If we set the diagonal of as the norm of the columns of , then is column stochastic, and if the input matrix is also column stochastic, then an exact solution requires the columns of to be stochastic as well, so that the columns of belong to the convex hull generated by the features in .
We now see how this considerations are important in practical applications.
1.2 Image Processing
One of the problem confronted by researchers in image processing is to decompose different images into common parts or features, both for identification purposes or for compression ones. For example, a common technique used in animation in order to contain the memory used is to not memorize into digital supports every pixel of each single frame, but to memorize only particular compressed or coded informations that lets a recorder to reproduce the film with little loss of quality.
In general, when confronted with a large set of images like the frames of a film, or a database of similar pictures, it can be convenient to memorize the common parts only one time, gaining space and also computational time for the recombining process. The problem is thus to find an efficient algorithm that automatically recognizes the common features and an intelligent way of storage of the informations.
Given a gray-scale image expressed as a matrix of pixels, with values in the real range , we can transform it into a real vector with as many coordinates as the pixels in the image. In particular, if , then we stack the columns of the matrix on top of one another, and obtain the vector defined as
Given a set of pictures of the same shape, we can now vectorize them and stack the corresponding vectors as the columns of our data matrix , and if we call the number of pixels of a single picture, becomes a nonnegative matrix in , so, after having fixed the number of common component we want to find, the NMF framework produces two matrices such that .
As already noticed, each column of is approximated by a linear combination of the columns of , that are nonnegative vectors of length . After having normalized by multiplication with a diagonal positive matrix (as discussed above), we can see its columns as images in the shape , so a generic column of , that is one of the original images, is now approximated as the superimposition of the pictures represented by some of the columns of .
Ideally, the images in are parts of the pictures in , like localized objects in the 2D space, so they’re usually sparse and disjoint images, that translates into sparse and nearly orthogonal vectors. In a famous experiment, Lee & Seung  processed a set of faces and the NMF automatically recognized their principal features like eyebrows, lips, noses, eyes, and so on, so that they were immediately human-recognizable. This example shows the importance of NMF as a decomposition tool for graphical entities.
As already said, the sparseness and the choice of are important factors. The sparseness is an index of the uniqueness of the solution, that is important on the side of interpretation of the output, since different solutions usually brings up set of pictures not human-recognizable as real objects and features. On the side of compression, we can see that the original pixels of are now coded into pixels in and coefficients in , so the compression is useful when the approximation is good with a low . On terms of images, it means that there are few components that span the whole set of pictures.
1.3 Transformations Issues
When we use NMF on a matrix
we usually expect the original images to have some predominant common features, so that the algorithm can find them with little noise. This may be true in the case of sets of static pictures, when calibrated and centered, but even in the case of facial recognition, there may be cases of misalignment, as already noticed by and many others. In general, the NMF suffers in this cases since it is not invariant under a vast set of transformations, for example shifts, rotations, symmetry, stretches and so on, in fact the common features must be in the exact same positions on the different pictures in order to be pinpointed.
This is a common problem faced in the animations programs, since, even if the subjects in a scene of a footage are the same, they constantly move on the screen, so their detection must follow some temporal scheme, and can’t be performed by a simple NMF.
Possible ways to deal with this problem are to change the data in one of the three matrices or . For example, if we add ta a transformed copy of each original picture for every transformation in a set we choose, then the common features get detected even if they’re deformed, but this increases the size of the problem by the square of the number of alterations used, that’s usually greater than the number of pixels in a single image. One possible solution is obviously to rise the parameter , but this leads to instability in the solution, as we already discussed.
A good idea seems instead to rise the quantity of data contained in the matrix , since we strife to maintain the graphical property of the columns of to represent the common features of the original images. In the next chapters we’ll define new notations and operators to deal with a matrix whose elements are capable to transmit more informations on pictures than simple real numbers.
In this document, our focus is on the problems related to the lack of translation invariance of NMF, so we’ll use shift permutations to modify the kind of elements contained in the matrix . First of all, we define an operator between matrices not necessarily real.
2.1 Diamond Operator
Given an element , represented by a couple , where is a real number, and is a permutation of indexes (that is, an element of the permutations group ), then it’s well defined its action on a real vector
The action of on makes it a linear operator, so it can be represented by a matrix, and in particular, since the action of each permutation is associated with a permutation matrix , it’s easy to see that
The algebra generated by the permutation group over the real field is denoted as , and its elements are finite sums of elements
As before, these elements have a natural action on , that is an extension of the action of , given by
so there exists an homomorphism of algebras that associates to each element of the algebra a real matrix, and later we’ll see how it behaves on a particular subgroup.
Let’s now suppose that is a matrix with entries in the above described algebra , and is a real matrix. We need an operator to apply the elements of to the columns of , so we define the diamond product :
Definition 2.1 (Diamond Product).
The diamond operator between a real matrix and a matrix is defined as
and returns a real matrix in .
In other words, the -th column of the diamond product is a linear combination of permutations of columns, with coefficients and permutations described by the elements of the ’s -th column.
Let’s also define the multiplication between two matrices with entries in the algebra of permutations. Remember that is an algebra, so sum and product are well defined, and the elements of can be viewed as well as matrices through the homomorphism , so the two operations correspond to the usual sum and composition of matrices.
Definition 2.2 (Diamond Product).
The diamond operator between two matrices and is defined as
and returns a matrix in .
This operation differs from the normal multiplication of matrices only because isn’t a commutative algebra, so we need to specify the order of the multiplication between the elements. The inverted order is necessary to partially maintain the associativity of the operation: given a real matrix , and two matrices with elements in the algebra, it’s easy to verify that
Ideally we need to invert the elements of and since is the first to act on the columns of , followed by .
One downside of this operation is that it doesn’t cope well with the normal matrix multiplication: given real matrices, and a matrix in the permutation algebra, then
Let’s now return to image transformations, and focus on a particular subgroup of the permutation algebra.
2.2 Shifts and Circulant Matrices
Given a gray-scale image , we’ve seen how to transform it into a vector . We want now to codify a shift on the image as a vectorial transformation: a shift of the original image by position on the horizontal axis and position on the vertical one will be encoded as a circular shift on of magnitude , that is, we produce a vector whose -th coordinate is the -th coordinate of .
If we call , we can denote as the cyclic subgroup of the permutation group whose elements shift cyclically all the indexes of vectors in by an integer constant. We’ll call the shift by position, where :
where the indexes are to be considered modulus .
The elements of are linear operators, so can be represented by matrices through the above mentioned homomorphism . In particular, the element is associated to the circulant matrix that has 1 on the first cyclic superdiagonal and 0 anywhere else, and , so that has 1 on the -th cyclic diagonal and zero otherwise.
In the next section, we’ll use elements of type to define a new problem with the same shape of a normal NMF, but on different domains, and since the shift is completely identified by the remainder class , we’ll refer to as the couple .
Now we reconsider the classic NMF, and widen the domain of the matrix . Our aim here is to find a new method to decompose pictures into common components, even when they’re shifted, so, like in the NMF, we stack the original images as columns of the matrix , and look for a matrix whose columns are the wanted common features, and a matrix with elements in , so that it can tell us both the intensity and the position of each component in into each original picture in .
In particular, we want to rewrite the NMF problem as
Problem 2.1 (PermNMF).
Given a matrix is in , we want to find a matrix in and a matrix in that minimize
The diamond operator is defined on elements of , but we restrict the entries of to elements in , so that a single image (column of ) is a linear combination of the images represented by the columns of , but shifted. We notice that expanding further the domain of usually leads to trivial and useless solutions; for example, if we let the elements of be in , that are linear nonnegative combinations of permutations in , then even with there’s a trivial solution that decomposes perfectly the matrix :
In other words, a linear combination of the translations of a single pixel can reconstruct any image, so it is an exact and completely useless solution. Moreover, expanding to the group usually leads to the dismembering of the images represented by the columns of , so we stick to work with this framework for this document.
An other particularity of this formulation is that, if we impose that each element of must be of the type , that is, we fix all the permutations to be the trivial identity, then the problem returns exactly the original NMF, and the diamond operator coincides with the normal matrix multiplication.
The PermNMF has the same structure of the normal NMF, so we can try to use similar solving algorithms. A characteristic we’d want from our solution is the sparsity of the columns, since they should represent isolated objects in the images, so the first algorithm considered is the MU update, since it is known to naturally produce sparse solutions. Unfortunately, the MU method efficiency, in the NMF case, comes from the approximation
but in our case, as already stated, there’s no associative property
For this reason, we resort to an ALS/PG setting.
The update of requires to solve a convex problem, so we can use some of the usual methods, like a modified Projected Gradient; this one is particularly good for this case, since we can’t transpose the expression in order to obtain the setting of the Active Sets algorithms.
For simplicity, we use the following PG algorithm, where we stop in case of low error or small step:
We’ll refer to this function from now on as
The computation for the gradient in the algorithm are developed in Appendix A, and it shows that
The operations performed in each cycle of the method have a computational cost of .
Let’s now focus on the update of , that requires to solve an optimization problem on the group . We start by solving a largely simplified problem.
3.1 Single Permutation NNLS
Let’s suppose to have two vectors in , and we want to find the best element of that minimizes
where the norm used is the euclidean one.
A natural assumption is that , otherwise every element gives the same value of . If we knew the optimal , then we could find without fail, because it becomes a simple Nonnegative Least Squares (NNLS) problem.
A simple solution consists into computing the optimal for every , and check which couple gives us the minimal error. We know that , so we can compute the error as a function of
The problem is thus equivalent to maximize , but we’re interested only in the positive case, so we focus on maximizing the scalar product , since if then for every , so .
By definition, is the vector shifted, so we can call the real nonnegative matrix that has all the shifted versions of as columns, and compute the maximal component of . Since is a circulant matrix, this operation costs
if performed with Fast Fourier Transformations, so this method is fast and gives us the correct solution.
From now on, we’ll use this algorithm with the syntax
Let’s now increment the number of permutations needed.
3.2 Multiple Permutation NNLS
Given now a vector , and a bunch of vectors we can now try to find the best elements that minimize the quantity
We’re thus looking for the best linear combination with positive coefficients of the shifted vectors that gives us the original vector . If we call the matrix with as columns, and the (column) vector of , then we can rewrite the problem in a compact way as
A way to solve this problem is using the precedent algorithm in an alternated fashion. In fact, if we fix , then it becomes a Singular Permutation NNLS problem on , and we know how to solve it exactly.
So we can solve the problem sequentially for each and repeat. The initial value of is usually given as an input parameter, but it can also be generated casually at the beginning of the algorithm.
From now on, we’ll use this algorithm with the syntax
Its computational cost is the number of iterations multiplied times the cost of The Single Permutation Problem, so it is considering as a constant. In particular cases, it may be useful to randomize the choice of the index , since it’s important not to impose a preference order on the components in .
3.3 Final Method
We can now return to the original problem
Like the normal NMF, it can be decomposed into smaller problems
that can be solved with the Multiple Permutation NNLS algorithm. If we put everything together, we obtain the final method
Every step of This ALS Update Method costs if we consider the number of iterations in the internal methods as constants. We will stop the updates when the convergence is too slow, when we loop on the same matrices, or when we reach a number of iteration too high.
3.4 Extension and Other Works
Given a set a pictures, now we’re able to perform a PermNMF and obtain a set of common features that can reconstruct the original data once combined through coefficients and permutations codified in . Given one of the images in , the algorithm tells us if it is present in the original images, but it doesn’t detect if it appears multiple times. One example of such instance may be a set of radar images, in which different objects intercepted by the wave signals have distinct shapes, but each one can appear multiple time in the same picture.
One possible solution is to perform an initial PermNMF with a parameter proportional to the effective number of distinct objects with multiplicity that can appear on a single image, discard the found components with low coefficients, and repeat the PermNMF on the output components with a low corresponding to the number of distinct shapes without multiplicity. Let’s call the first larger parameter, and the set of pictures to analyze. We obtain
where , and , , so the final decomposition will be again a real matrix with components, and a matrix . This last matrix is able to tell, for each component, even if there are multiple instances in every original image.
The computational cost of such method (for each cycle, till convergence) is
that, under the assumption , is equivalent to , meaning that the second step has a negligible computational cost compared to the first. If is still on the order of magnitude of , the asymptotic cost doesn’t change, but if that’s not the case, it is better to look for other ways.
On this topic, Potluru, Plis and Calhourn in  offer an algorithm that uses Fast Fourier Transformations and circulant matrices in order to compute and codify permutations of the components, called ssiNMF (sparse-shift invariant NMF). As in the PermNMF, the basic idea is to find components and a set of permutations that could reconstruct the original images, but the ssiNMF sets as target the permutations in the group , corresponding through with all the circulant nonnegative matrices, so that all the operations can be performed through FFTs. Thanks to this, their algorithm is able to directly construct an approximation
Eggert, Wersing and Korner in  took a more general approach to the problem: as we set a subgroup of , they chose a general set of transformations of the plane, seen as operators on the columns of , and multiplied the number of parameter of by the cardinality of the chosen set, so that for each transformation of the components there would be coefficients in stating their intensity in the original images.
Both the approaches suffer by the presence of the trivial and exact solution described in section 2.3: a single pixel can generate any image if we allow too many transformations of the space. They propose to perform a common modification on the NMF framework, that is adding a penalty factor to ensure the sparseness of the output, since the presence of a single pixel in the component output corresponds to a lot of positive coefficients in , and it leads to the presence of an additional parameter to set manually or through validations techniques.
An other characteristic of both the algorithm is the rise in memory used and asymptotic computational cost by at least a factor on par with the number of pixels on a single image, leading to a cost by iteration at least of . When compared with the PermNMF algorithm, we see that they’re comparable when , meaning that a component have to appear in the original image on average times.
In these experiments, we use the PermNMF algorithm seen in the previous chapter, with the initial parameters and generated randomly, and the variable set to 10 in both the and the methods.
In the first experiment (Figure 1) we use 2 simple shapes (a square and a cross) of 9 pixel that move into a frame of dimensions 20x20, and add a casual error of mean 0.15 (where each pixel has an intensity between 0 and 1). In this case the algorithm manages to find the right components after less than 10 repetitions on average. The images shown on the bottom row are the column of , and they’re distinguishable as a cross and a square, with little noise given by the imperfections on the original images.
In the second experiment, we generate 20 images of shape 30x30 from three simple figures (a plane, a tank and a ship), with a nois of mean 0.15. Each image can include up to two copies of the same figure, so we need to perform a first PermNMF with , and then a second time with to extract the original ones. The first application of the algorithm is slowed down by the presence of the same shapes multiple times in the images, but the second application is real fast. As said, we managed to extract first the common features with their multiplicity, and then the actual features. Multiplying the two matrices we obtained in the two steps of the algorithm, we can deduce the actual position with multiplicity of the shape found in all the 20 original images.
5 Future Works
The PermNMF has not been throughly studied and analyzed. First of all, it lacks a convergence result, both because the usual arguments used for the ALS algorithms vastly use the fact that the two subproblems in the classical NMF are convex, and because we switched the framework to non-continuous spaces such as , where it is still not even well defined a canonical concept of "local minimum" (the usual topological embedding of this space in gives a notion of stationary points that doesn’t cope well with the nature of permutations).
On the point of view of the PermNMF problem, there’s a lot to say, for example, on whether there exists an exact algorithm, or if there are bounds on the minimum , or even if the solution is unique (up to trivial transformations). In , Gillis find a preprocessing for the input data that gives a more well-posed problem then the normal NMF, so such a transformation could be beneficial even to the PermNMF. In , the authors found precise conditions for under which there exists a polynomial time algorithm for the exact NMF problem, and stated that in general the approximation problem is NP-hard, so it’s highly possible that even the PermNMF problem is a NP-hard problem, and that a the polyomial time algorithm could be adapted for this case.
On the side of the algorithm itself, it’s possible that a MU (Multiplicative Update) approach on , even if expensive, could retain its descend property, so it can become a substitute or an aid for the PG method. On both the update of and , it is still possible to apply a CD (Coordinate Descend) method, even if it also lost most of his efficiency due to the bad behavior of the diamond operator. Both this methods, MU and CD, are also recommended for the generation of sparse solutions, a feature we’d like to obtain. On the Multiple PermNNLS algorithm, moreover, it’s also possible to consider an active-set like method to choose preemptively which element to update in every cycle, in order to make the error drop faster.
Eventually, we studied the problem when the elements of are restricted to , but it’s possible also to consider other subgroups and subalgebras of in order to encode different transformations of the plan, or just to make the NMF invariant with respect to particular linear operators.
Appendix A Computation of gradient
Let’s compute the gradients needed.
In the following steps, we consider the general element of as a (circulant) matrix, using implicitly the homomorphism .
If we denote the matrix as the couple , then its transpose is represented by the couple . Let’s call the matrix with the same dimension of and , so we have
We can continue the computation as
So we can write in a compact form the gradient
Sanjeev Arora, Rong Ge, Ravindran Kannan, and Ankur Moitra.
Computing a nonnegative matrix factorization – provably.
Proceedings of the Forty-fourth Annual ACM Symposium on Theory of Computing, STOC ’12, pages 145–162, New York, NY, USA, 2012. ACM.
-  J. Eggert, H. Wersing, and E. Korner. Transformation-invariant representation and nmf. In Neural Networks, 2004. Proceedings. 2004 IEEE International Joint Conference on, volume 4, pages 2535–2539 vol.4, July 2004.
-  Brendan Frey and Nebojsa Jojic. Transformation-invariant clustering and dimensionality reduction using em. IEEE Trans. Pattern Analysis and Machine Intelligence, 2000:1–17, 2000.
-  Nicolas Gillis. Sparse and unique nonnegative matrix factorization through data preprocessing. J. Mach. Learn. Res., 13(1):3349–3386, November 2012.
-  Daniel D Lee and H Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401(6755):788–791, 1999.
-  Vamsi K Potluru, Sergey M Plis, and Vince D Calhoun. Sparse shift-invariant nmf. In Image Analysis and Interpretation, 2008. SSIAI 2008. IEEE Southwest Symposium on, pages 69–72. IEEE, 2008.