FRADM
None
view repo
We address the problem of efficient sparse fixedrank (SFR) matrix decomposition, i.e., splitting a corrupted matrix M into an uncorrupted matrix L of rank r and a sparse matrix of outliers S. Fixedrank constraints are usually imposed by the physical restrictions of the system under study. Here we propose a method to perform accurate and very efficient SFR decomposition that is more suitable for largescale problems than existing approaches. Our method is a grateful combination of geometrical and algebraical techniques, which avoids the bottleneck caused by the Truncated SVD (TSVD). Instead, a polar factorization is used to exploit the manifold structure of fixedrank problems as the product of two Stiefel and an SPD manifold, leading to a better convergence and stability. Then, closedform projectors help to speed up each iteration of the method. We introduce a novel and fast projector for the SPD manifold and a proof of its validity. Further acceleration is achieved using a Nystrom scheme. Extensive experiments with synthetic and real data in the context of robust photometric stereo and spectral clustering show that our proposals outperform the state of the art.
READ FULL TEXT VIEW PDF
In this work, we address the problem of outlier detection for robust mot...
read it
Matrix learning is at the core of many machine learning problems. To
enc...
read it
Spectral clustering is a very important and classic graph clustering met...
read it
Kmeans  and the celebrated Lloyd algorithm  is more than the cluste...
read it
In this article, we derive a Bayesian model to learning the sparse and l...
read it
This paper studies clustering for possibly high dimensional data (e.g.
i...
read it
Treestructured data usually contain both topological and geometrical
in...
read it
None
Systems with fixedrank constraints exist in many applications within the fields of computer vision, machine learning and signal processing. Some examples are: photometric stereo, where depth is estimated from a still camera that acquires images of an object under different illumination conditions, leading to a rank constraint; motion estimation, where the type of motion of the objects defines a rank.
This paper addresses the problem of efficient sparse fixedrank (SFR) matrix decomposition, i.e.: given a matrix affected by outliers, this is, gross noise of unknown magnitude, we aim to recover an uncorrupted matrix and a sparse matrix such that and , with known beforehand, as defined in (1),
(1) 
SFR matrix recovery is intimately related to the sparse lowrank (SLR) recovery problem (2
), for which algorithms such as Robust Principal Component Analysis (RPCA)
[3] and Principal Component Pursuit (PCP) [30] are well known due to their extraordinary capabilities to solve it and their application to a wide range of problems.(2) 
Robust SFR recovery might seem a simpler case of SLR decomposition, or even a straightforward derivation. However, SFR recovery is a hard problem that involves a highly nonconvex constraint due to the rank imposition. This factor is not present in the SLR decomposition problem due to the nuclear norm relaxation. Therefore, a careful design is needed in order to produce a stable SFR decomposition method with a good convergence rate.
In addition to the convergence speed, achieving efficient and scalable SFR decompositions requires algorithms with very low computational complexity per iteration. The main bottleneck of these algorithms is the enforcement of the correct rank or its minimization, a step that usually requires the use of a TSVD or an SVD per iteration, which complexity is for a matrix of rank . How to reduce this bottleneck is a line of research that has been recently targeted by several works such as [9][23][24], showing interesting ideas leading to algorithms with quadratic and linear complexities with respect to the input matrix size. The key lessons to learn from these works are two: (i) the factorization of largescale problems into products of small size matrices [16]; and (ii) the use of a subsampled version of the input matrix to produce fast and accurate approximations of the solution [24].
Our work has been influenced by these concepts and several ideas drawn from stateoftheart differential geometry techniques. We have experimented with the mentioned concepts and improved upon them in order to create an efficient and precise SFR decomposition algorithm suitable for large scale problems. In this respect we present the following contributions: (i) an optimization method, named FRADM^{1}^{1}1Code is available at https://github.com/germanRos/FRADM (FixedRank Alternating Direction Method), that solves SFR problems following an ADM scheme to minimize an Augmented Lagrangian cost function; (ii) a novel procedure to impose fixedrank constraints through a very efficient polar factorization, named FixedRankOptStep, which is superior in convergence, stability and speed than the bilinear counterparts used by stateoftheart methods and; (iii) the use of a simple projector to impose SPD constraints efficiently along with a novel proof of its validity.
We show that our method, based on the FixedRankOptStep procedure, outperforms in time, accuracy and region of applicability current stateoftheart methods. Furthermore, we also show that our proposal FRADM can benefit from Nystrom’s method [25] to improve its computational efficiency while maintaining a good level of accuracy. These results are supported by thorough experimentation in synthetic and real cases, as detailed in Sec. 5.
Capital letters, such as
represent matrices, while vectors are written in lowercase.
stands for the matrix transpose, for its pseudoinverse and is the matrix trace operator. stands for theth largest singular value of a given matrix. The indexation of the
th row and the th column is defined as . Matrix subblocks of are referred to as to index from row to and column to . is the Frobenius norm and are the entrywise norm and the matrix nuclear norm, respectively. and are the square and the rectangular identity matrices. , is the Stiefel manifold of matrices with . and stand for the Symmetric (Semi)Positive Definite matrices, respectively. is the fixedrank manifold of matrices with and is the set of fullrank matrices. stands for the Orthogonal group, but be careful, since is also used to describe the complexity of algorithms in bigO notation. We also make use of some proximity operators and projectors defined as: , the symmetric part of M. for the standard softthresholding (promotes sparsity); for the projector onto the Stiefel manifold, and for the projector onto the SPD manifold (these are defined in Sec. 4.1).The Accelerated Proximal Gradient (APG) [13] will serve as our starting point within the plethora of methods present in the literature. This method, although it is not the first one proposed to solve the RPCA problem (see for instance FISTA [2]), is an appealing combination of concepts. It approximates the gradients of a given cost function to simplify its optimization and improves its convergence. It also includes Nesterov updates [20] and the critical continuation scheme [13], which all together lead to a method with sublinear convergence , where is the number of iterations. Its computational complexity per iteration is for matrices. Afterwards, authors of [12] proposed the Augmented Lagrangian Multiplier method (ALM) in two flavours. First they present the exact ALM (eALM), which uses an Alternating Direction Method (ADM) to minimize an Augmented Lagrangian function in a traditional and exact way. Then, an inexact version is also proposed (iALM), which approximates the original algorithm to reduce the number of times the SVD is used. The convergence rate of eALM depends on the update of , the penalty parameter of the Augmented Lagrangian. When the sequence grows geometrically following the continuation principle, eALM is proven to converge Qlinearly . For iALM, there is not proof of convergence, but it is supposed to be Qlinear too. Both methods have computational complexity of per iteration.
Recently, ALM was extended in [14], which included a factorization technique along with a TSVD from the PROPACK suite [11] to achieve a complexity of per iteration. The bottleneck caused by the TSVD has also been addressed via random projections, leading to the efficient RandomizedTSVD (RTSVD) [9]. However, although being more efficient than the regular TSVD, results are considerably less accurate. The idea of including a factorization of the data was then improved by LMAFIT [23], which uses a bilinear factorization to produce two skinny matrices and , such that , to speed up the process. A similar concept was used in the Active Subspace method (AS) [16], but in this case the bilinear factorization is given by and , such that . This formulation turns out to be very useful when , leading to a complexity per iteration of . Unfortunately, is an upper bound for the actual rank of and needs to be given by the user. This is not suitable for SLR scenarios, but fits perfectly in the SFR framework. Another point to highlight about LMAFIT and AS is the utilization of closedform projectors to impose constraints like orthogonality, lowrank and sparsity. This algebraical way of optimizing functions differs from the geometrical counterparts in the literature on manifold optimization (see [1][18]). Substituting all the required machinery to perform differential geometry (e.g., retractions, lift maps, etc.) by projectors seems a good idea from the point of view of efficiency. However, this method is not absent of problems. The factorization in AS is highly nonconvex, an issue that influences the number of iterations required for the convergence of the method, which is notably higher than eALM, despite having the same theoretical convergence rate .
One of the contributions of our work is to improve the convergence of fixedrank projection methods. To this end we employ a polar decomposition as in [18]. This polar decomposition offers us the possibility of exploiting the manifold structure of fixedrank problems as the product of two Stiefel and an SPD manifold. . However, we deviate from [18] to propose more efficient expressions that make use of projectors to speed up the process, giving rise to a better convergence. We also consider worth highlighting a key tool described in the recent work [24]. There, the authors follow a strategy that resembles the one described in [16], but they add a subsampling step based on the Nystrom’s method [25] that leads to a linear complexity per iteration. We borrow this idea to further speed up our optimization.
We propose the resolution of the nonconvex program in (3) as a direct way to perform the sparse and fixedrank decomposition—note that (3) is equivalent to the program (1) defined in Sec.1.
(3) 
The optimization of (3) is carried out over an Augmented Lagrangian function, leading to (4). stands for the Lagrange multiplier and represents the fixedrank manifold of rank .

(4) 
To efficiently solve (4) we utilize an ADM scheme [12] endowed with a continuation step, as presented in Algorithm 1. The update of the fixedrank matrix is obtained via the FixedRankOptStep algorithm that implements the proposed polar factorization. For the sparse matrix , the standard softthresholding is used. Notice that is updated following a geometric series (Alg.1.#9) in order to achieve a Qlinear converge rate [12]. Despite having the same asymptotic convergence order as LMAFIT, AS and ROSL, our method FRADM takes less iterations to converge, due to the accuracy of the novel FixedRankOptStep. This is especially important for the cases where the magnitude of the entries of are similar to those of , a challenging situation that other stateoftheart methods fail to address correctly. We provide empirical validation for this claim in Sec. 5.
An adapted version of FRADM, referred as FRNys, is also provided. FRNys exploits Nystrom’s subsampling method ([25][24]), to further speed up computations. This method is presented in Algorithm 2 and follows the recipe given in [24].
Nystrom’s scheme proceeds by randomly shuffling the rows of producing . Then, the top and the left blocks of are processed separately by using FRADM (Alg.1). Notice that these blocks are chosen of size and , where has to be a number larger than the expected matrix rank. In our case . Finally the independently recovered matrices and are combined to produce and .
Imposing rank constraints requires an efficient way of computing the projection of an arbitrary matrix with arbitrary rank onto the fixedrank manifold . A simple solution is provided by the EckartÃ¢â‚¬â€œYoung theorem [5], which shows that the optimization problem (5):
(5) 
is solved by the truncated SVD (TSVD) of . Despite the success of the TSVD as a tool for producing lowrank approximations, and the many available improvements, as for instance the usage of random projections[9], some problems require the computation of many TSVDs (typically one per iteration) of very large matrices. Thus an efficient alternative to the usual TSVD algorithm is required.
In this section we propose the method FixedRankOptStep (Algorithm 3), which computes a fast approximate solution to the projection of a matrix onto the fixedrank manifold, like the given by the TSVD but much faster. Additionally, in the Appendix we also propose the method FixedRankOptFull (Algorithm 4) that can be seen as a series of iterations of the FixedRankOptStep algorithm, providing a solution with a prescribed accuracy and with Qlinear convergence rate to the minimization problem (5) living on . The FixedRankOptStep algorithm is suitable for largescale problems where many TSVDs of large matrices are required, and an approximate solution is faster and enough for convergence, as we will show later.
Following [18], we use a polar factorization on suggested by the TSVD. Given a matrix of rank , its TSVD factorization is
(6) 
where , and . Then, a transformation
(7) 
where , does not change , and allows to write it as
(8) 
where now , , and . Thus, the fixedrank manifold can be seen as the quotient manifold . From this, we reformulate (5) in as the solution of (9).
(9) 
The FixedRankOptStep algorithm performs a single step of an alternating directions minimization (ADM) on each of the submanifolds , and (Algorithm 3).
The minimization subproblems on Stiefel manifolds involving and in Algorithm 3, are not the standard Stiefel Procrustes Problem [6]. Here, the Stiefel matrix is leftmultiplying instead of rightmultiplying, as usually, which allows to provide a fast closedform solution by using the Orthogonal Procrustes Problem (OPP) [22], as shown in (10):
(10) 
where denotes the projector onto the Stiefel Manifold. This can be efficiently computed through a skinny SVD as . Alternatively, if (maximal rank, as we shall assume in the following), it can be computed as . This shows that always exists and it is unique. A similar result holds for the minimization of by simply transposing (10).
The minimization subproblem on the SPD manifold is more challenging. The reason is that, although convex, the SPD manifold is an open manifold and therefore the existence of a global minimum is not guaranteed. Its closure is the SPSD manifold, and there the existence of a solution is neither guaranteed. However, we shall see that in our case there exists a minimun in . Let us analyse this by first introducing a novel projector onto the SPD manifold. To this end we consider the SPD Procrustes Problem [26] (11):

(11) 
where the projector is simply given by . In general, the solution of the SPD Procrustes Problem requires solving a Lyapunov equation [26], but in our case is simpler since and are orthogonal. Although in general there is not guarantee that is positive definite, we can assure it for our formulation, see the Appendix.
Since the optimization problem (3) is highly nonconvex, a global convergence theorem as in eALM [12] cannot be given. However, a weak covergence result similar to iALM or that of LMAFIT [23] (where there is no nuclear norm minimization) can be given. For that purpose, let us state the firstorder optimality conditions for the constrained minimization problem (3):
(12)  
where , and is a Lagrange multiplier. Then we can prove that:
Proof: Using Algorithm 1, and given a projector , we have:

FRADM, and its Nystrom accelerated variant FRNys, are compared here against the selected methods of the state of the art, i.e.: Accelerated Proximal Gradients (APG) [13], inexact Augmented Lagrangian Multiplier (iALM), exact Augmented Lagrangian Multiplier (eALM) [12], Active Subspace method (AS) [16], LowRank Matrix Fitting (LMAFIT) [23] and Robust Orthonormal Subspace Learning (ROSL) [24]. These methods are good representatives of the evolution of SLR and SFR solutions, ranging from the fundamental proximal gradients of APG to sophisticated factorizations included in AS and ROSL, via ADMM optimization [12]. We also included a version of iALM that makes use of the Randomized TSVD (RTSVD) [9] in order to show the benefits of our approach against simple randomization. It is critical for the correct understanding of our experiment to clarify that we have split the previous methods up into two categories, i.e., SLR and SFR techniques. APG, iALM and eALM represent SLR methods, i.e., the rank is not known a priori; while RTSVD, AS, LMAFIT, ROSL and FRADM represent SFR solutions, i.e., a correct initialization of the rank is provided, since the specific application allows it. This assumption holds for the entire section. Experiments are conducted over synthetic and real data to show the capabilities of our technique in computer vision problems. All the algorithms have been configured according to the suggestions of their respective authors. The experiments were run on a desktop PC at 3.2GHz and 64GB of RAM.
We test the recovery accuracy and time performance of the methods with matrices of different dimensions and ranks. To this end, we generate fullrank matrices and
from a Gaussian distribution
, such that and . A sparse matrix representing outliers is created with a given percentage of its entries being nonzero and magnitudes in the range . Then, the final corrupted matrix is . We deliberately forced the sparse entries to have a magnitude similar to the one of the expected lowrank matrix. The reason for this is that usually the experiments presented in the literature impose a good differentiation between the magnitude of the entries of and , making the recovery problem almost trivial. Here, we remove that simplification, allowing for similar magnitudes of the corrupted entries, which makes the problem more interesting. We will show that, with this challenging setup, the performance of many stateoftheart methods dramatically decreases, while our approach maintains a good recovery accuracy.Our first test measures the recovery capabilities of the different methods under study when subjected to similar magnitudes of the entries of and . To this end we create corrupted matrices of increasing rank and an increasing fraction of outliers. The result of this experiment is shown in Figure 2 in form of phase transition diagrams, with rank fractions represented in the xaxis and outlier fraction in the yaxis. Colors represent the recovery (inverse) probability of each case, i.e., the lower error (cold colors, i.e. blueish) the better. From this plot, it can be seen that these conditions are very challenging for all the algorithms.
APG, eALM and iALM, making use of an SVD, end up with a very narrow recovery region (in blue). RTSVD gets a narrow recovery region due to accuracy problems (see also Fig. 1 for further information). Notice that AS is not even able to converge beyond a of rank due to the strong nonconvexity induced by its bilinear factorization. LMAFIT shows a rather acceptable recovery region, while ROSL clearly suffers in obtaining a correct recovery for this sort of data. In our analysis, ROSL performs well when the magnitude of (the noisy entries) are several magnitudes bigger than those of . However, in other cases the recoverability of ROSL dramatically decay. Our proposal, FRADM, presents the best recovery for a wider region even in this challenging setup. This characteristic is critical for real applications where outliers might be either very large or very subtle.
We also evaluated one of the most critical aspects of these methods, i.e., the accuracy of a given method at providing a good lowrank approximations of a matrix . Stateoftheart approaches have gained in efficiency by replacing the SVD for a more convenient fixedrank projection, as in the case of RTSVD, AS, LMAFIT, ROSL and our proposal FRADM. However, as shown in Figure 1, different projection strategies lead to different convergence rates and speeds. In this way, when compared against an exact TSVD, the polar decomposition used by FRADM turns out to be superior to all its competitors, as derived from the reduced number of iterations required to achieve a relative error of . We would like to highlight that our approach even presents a better convergence behaviour than the wellknown RTSVD, which is considered one of the fastest methods for lowrank projection. Later, we will show that FRADM not only has a better convergence, but is also faster and more accurate.
Our second experiment uses matrices of increasing sizes (), ranging from to , while keeping the rank fixed, and the entries magnitudes as defined above. 10 repetitions are considered per each size. In this case the methods under evaluation are the APG, iALM, eALM, RTSVD, AS, LMAFIT, ROSL, ROSL+ and our proposal FRADM, along with its equivalent accelerated version, FRNys. We have accelerated FRADM to present a counterpart to ROSL+[24]. In this way we can offer a fair comparison with our proposal and show that our method remains superior after the Nystrom’s speedup. Results of this test are shown in Table I, considering the recovery error for both matrices and given by and , where and are the optimal matrices. We also consider computational time in seconds and the number of iterations used by each method. FRADM is the method with the best tradeoff of high recovery accuracy and low computational time (the fastest of the nonaccelerated methods). The efficiency of methods such as LMAFIT and ROSL considerably decreased due to their difficulties to face small sparse entries. APG, iALM and eALM also find troubles searching for the appropriate rank in this challenging conditions. For the case of the RTSVD, its accuracy is lower than desired, and due to its lack of accuracy requires too many iterations to converge.
For the accelerated methods, FRNys has proven to be the fastest and the most accurate in all synthetic tests, despite the accuracy degradation provoked by the matrix sampling. The benefits of applying Nystrom’s acceleration are clear, specially for big matrices, as in the case, where the total time is reduced in two orders of magnitude. However, as we show in the next experiment, this acceleration is not convenient for problems with large matrix ranks.
In this third experiment methods performance is tested against matrices of increasing dimensions and rank. Matrices are created as described above, but their rank is established to be , where is the matrix size. Results are summarized in Table II. The first thing to notice is that the time of Nystromaccelerated methods is bigger than their unaccelerated counterparts. This is due to the high rank of the problem and that the matrices resulting from the Nystrom’s sampling technique are of sizes and , with big enough (usually ). This leads to two matrices that are almost of the size of the original one, making the use of Nystrom counterproductive. Regarding the unaccelerated methods, FRADM performs almost twice faster than the second best approach, AS. In terms of recovery accuracy, all the methods present similar results, except for small matrices, where ROSL fails due to its sensibility to initialization parameters.
SLR  SFR No Accel.  SFR Accel.  
APG  iALM  eALM  RTSVD  AS  LMAFIT  ROSL  FRADM  ROSL+  FRNys  
Err. L  6.2e7  6.3e10  3.3e9  1.0e7  7.6e9  1.3e4  5.0e10  2.21e10  1.7e2  6.8e10  
Err. S  8.4e5  1.1e7  5.8e7  8.7e6  3.8e7  9.1e3  7.5e8  3.6e8  1.2e+0  4.8e8  
500  iters  140  33  9  96  120  40  93  28  200  68 
time  11.74  1.25  2.29  0.90  0.59  0.11  0.89  0.11  0.67  0.17  
Err. L  4.7e7  1.1e9  2.5e10  1.8e7  1.4e9  1.8e7  1.3e9  5.0e10  1.2e4  6.2e10  
Err. S  8.7e5  5.1e7  6.4e8  1.6e5  4.3e7  1.2e5  6.3e7  2.0e7  8.4e3  8.7e7  
1K  iters  142  34  10  95  133  65  98  28  200  65 
time  61.75  5.87  11.04  1.57  2.39  0.75  4.27  0.46  1.29  0.22  
Err. L  3.7e7  8.1e10  2.3e10  4.4e7  1.3e9  2.1e8  1.1e9  3.0e10  5.2e8  3.1e10  
Err. S  9.4e5  4.5e7  7.8e8  3.4e5  6.3e7  4.6e7  6.9e7  1.8e7  3.6e6  7.6e8  
2K  iters  144  35  10  92  131  300  98  29  200  66 
time  396.4  20.34  50.02  5.14  9.64  12.45  17.79  1.57  1.98  0.31  
Err. L  2.6e7  4.8e10  2.5e10  7.3e7  9.4e10  1.3e8  6.8e10  2.7e10  4.7e9  3.3e10  
Err. S  9.3e5  4.2e7  1.0e7  5.9e5  6.3e7  2.9e7  6.0e7  2.2e7  3.2e7  2.2e8  
4K  iters  147  36  10  91  135  300  99  29  140  62 
time  3002  112  328  18  50.05  56.93  67.63  7.52  3.33  0.65  
Err. L  1.9e7  5.4e10  2.7e10  1.4e6  6.7e10  9.4e9  4.9e10  4.5e10  5.4e9  5.0e10  
Err. S  9.3e5  6.5e7  1.6e7  1.3e4  6.8e7  2.0e7  6.2e7  7.6e7  3.7e7  2.8e8  
8K  iters  150  36  10  88  138  300  100  28  139  61 
time  22415  517  2214  63.6  243  202  261  27.42  6.90  1.24 
SLR  SFR No Accel.  SFR Accel.  
APG  iALM  eALM  RTSVD  AS  LMAFIT  ROSL  FRADM  ROSL+  FRNys  
Err. L  1.3e5  1.7e4  1.1e6  2.3e7  1.5e8  1.8e4  1.2e2  9.0e9  1.7e4  2.7e10  
500  Err. S  1.3e4  1.7e3  1.1e5  4.1e5  1.1e7  2.8e2  1.2e1  1.1e7  2.5e2  4.1e8 
r=50  iters  175  37  11  88  85  48  138  37  200  66 
time  13.77  3.61  15.44  1.70  0.75  0.42  6.43  0.34  8.99  0.97  
Err. L  2.4e6  6.9e9  4.4e7  6.7e7  1.6e8  2.8e8  1.2e8  9.9e9  4.3e8  8.1e10  
1K  Err. S  3.4e5  1.3e7  6.2e6  1.6e4  1.7e7  2.2e6  8.8e8  1.8e7  9.4e6  1.8e7 
r=100  iters  174  37  11  82  79  300  69  36  200  62 
time  68.88  14.36  61.74  7.63  3.18  7.43  33.13  1.54  69.75  3.61  
Err. L  6.8e7  6.3e9  2.0e7  1.7e6  1.1e8  1.9e8  1.1e8  1.1e8  4.5e8  1.6e9  
2K  Err. S  1.3e5  1.7e7  4.0e6  5.9e4  1.6e7  2.0e6  1.2e7  3.0e7  1.4e5  4.8e7 
r=200  iters  175  37  11  77  85  300  67  35  200  60 
time  461  80.49  332  32.74  17.32  30.77  325  7.53  653.13  16.93  
Err. L  3.5e7  7.5e9  1.5e7  5.0e6  1.2e8  1.3e8  1.0e8  6.8e9  4.7e8  5.3e9  
4K  Err. S  1.0e5  2.8e7  4.4e6  2.4e3  2.7e7  1.9e6  1.8e7  2.6e7  2.0e5  2.3e6 
r=400  iters  175  37  10  71  89  300  66  36  200  57 
time  3453  529  2106  183.5  107  162  3586  43  7008  96.57  
Err. L  5.9e7  5.0e9  4.3e9  1.5e5  6.3e9  8.6e9  7.1e8  9.5e9  4.8e8  1.4e8  
8K  Err. S  3.7e4  5.5e6  3.4e6  1.0e2  6.6e6  1.8e6  1.2e5  1e5  3.0e5  8.6e6 
r=800  iters  143  35  10  65  130  300  7107  21  200  54 
time  23130  2651  6394  1382  1075  1035  97397  166  91242  564 
We have chosen photometric stereo [27] (PS) as our first example of fixedrank problem. PS consists in estimating the normal map and depth of a still scene from several 2D images grabbed from the same position but under different light directions. The Lambertian reflectance model [27] is assumed, such that the light directions , the matrix of normals (unknowns) , and the matrix of pixel intensities are related via , where represents the albedo. The objective of recovering the normal map can be achieved by a LeastSquares (LS) method, but the quality of such a solution would suffer in the presence of outliers. Instead, robust decompositions can be used to get ride of outliers, as proposed in [28]. Since is a product of two rank3 matrices, in ideal conditions its rank is at most 3. We make use of this rank property to recover an uncorrupted version of that leads to a better estimation of the map and consequently of the depth map.
In our tests we use a dataset of objects viewed under 20 different illuminations, provided in [29]. From such images, we recover an uncorrupted version of the intensities . Then we run the Photometric Stereo Toolbox [29] to recover normal maps, depth maps, 3D models and some statistics. Table III shows the error in the normal maps after the recovery process with different methods. Here, we consider the reconstruction error, i.e., the normal map is rerendered into a shading image and then compared with the captured images. From the resulting error map several statistics are computed (RMS, mean and maximum error). The classical LS approach is taken as a reference of nonrobust approaches. As robust methods, APG, eALM and iALM, AS, LMAFIT, ROSL and FRADM are considered. RTSVD has not been considered due to its observed reduced accuracy. Nystrom accelerated versions are excluded due to the small size of the observation matrices, a constraint that prevents speedups.
The comparison shows that AS, ROSL and FRADM are the most accurate methods, producing estimations of the normal map with reconstruction errors below . The remaining methods are far from the accuracy offered by these fixedrank techniques, producing high residuals. Although AS consistently presents a lower error in the majority of the cases, the error differences below are of no impact for the application. This is shown in the error maps of Fig. 3(a). However, computational time is a critical factor for this problem, where, FRADM is one order of magnitude faster than ROSL and two orders faster than AS.
This figure displays the error maps of the considered approaches. As expected, LS leads to high errors due to outliers. APG, iALM and eALM improve LS results, but since they do not use the rank constraint recovered matrices have an erroneous lowrank. Fixedrank techniques, such as AS, ROSL and FRADM achieve very low residuals, making the error maps black. The recovered normal maps after the application of the FRADM technique are shown in Fig. 3(b) along with the 3D reconstruction of the objects. It can be concluded that SFR techniques, can drastically benefit problems like photometric stereo and FRADM stands as the fastest alternative while offering a very high accuracy.
LS  APG  iALM  eALM  AS  ROSL  LMAFIT  FRADM  

RMS  1.4e2  3.7e3  3.9e3  3.9e3  1.2e12  1.6e11  2.3e2  1.5e11  
Frog  Mean Err.  1.1e2  2.7e3  2.7e3  2.7e3  1.2e12  1.4e11  7.9e3  1.3e11 
Max Err.  1.6e1  2.2e2  2.1e2  2.1e2  1.8e12  4.8e11  2.1e1  4.7e11  
Time(s)  x  2.3e+2  1.4e+2  5.6e+2  3.1e+1  4.0e+1  1.4e+2  7.1e+0  
RMS  1.4e2  2.7e3  2.5e3  2.5e3  4.3e14  2.7e11  9.6e3  2.5e11  
Cat  Mean Err.  9.3e3  1.9e3  1.8e3  1.8e3  4.1e14  2.3e11  3.9e3  2.2e11 
Max Err.  2.2e1  1.8e2  1.4e2  1.4e2  6.4e14  6.6e11  1.4e1  6.7e11  
Time(s)  x  1.8e+2  1.1e+2  4.3e+2  2.4e+1  3.0e+1  1.1e+2  5.9e+0  
RMS  1.5e2  2.9e3  2.8e3  2.8e3  6.0e13  2.6e11  1.4e2  2.6e11  
Hippo  Mean Err.  9.8e3  1.6e3  1.5e3  1.5e3  5.7e13  2.4e11  6.4e3  2.3e11 
Max Err.  1.9e1  2.3e2  1.9e2  1.9e2  9.8e13  8.1e11  1.8e1  8.4e11  
Time(s)  x  1.9e+2  1.2e+2  4.7e+2  2.6e+1  3.2e+1  1.2e+2  6.0e+0  
RMS  1.4e2  4.0e3  3.9e3  3.6e3  3.8e12  1.8e11  1.8e2  1.5e11  
Lizard  Mean Err.  1.2e2  3.1e3  3.0e3  2.8e3  3.5e12  1.6e11  6.2e3  1.3e11 
Max Err.  1.7e1  3.6e2  2.7e2  2.7e2  1.2e11  5.5e11  2.2e1  4.4e11  
Time(s)  x  2.8e+2  1.6e+2  7.8e+2  3.7e+1  4.3e+1  1.6e+2  8.9e+0  
RMS  1.0e2  2.7e3  2.5e3  2.5e3  1.4e11  1.9e14  1.5e2  6.8e11  
Pig  Mean Err.  7.9e3  2.2e3  2.1e3  2.0e3  1.4e11  1.5e14  5.1e3  5.5e11 
Max Err  2.1e1  1.2e2  1.5e2  1.4e2  2.7e11  8.7e14  2.2e1  3.1e10  
Time(s)  x  2.3e+2  1.4e+2  5.2e+2  3.2e+1  3.7e+1  1.5e+2  7.7e+0  
RMS  4.3e2  1.1e2  9.1e3  9.9e3  8.8e13  2.8e13  2.7e2  1.3e13  
Scholar  Mean Err.  3.3e2  1.0e2  8.4e3  9.2e3  7.9e13  2.2e13  1.5e2  1.0e13 
Max Err.  3.3e1  3.3e2  2.2e2  2.4e2  1.9e12  1.3e12  2.4e1  6.0e13  
Time(s)  x  5.0e+2  3.0e+2  1.3e+3  6.5e+1  8.0e+1  3.1e+2  1.5e+1 
We address clustering as a fixedrank optimization problem with a known number of clusters represented by the matrix rank, where such a rank can become very high. Here, SFR methods can be easily added to the pipeline of Spectral Clustering approaches (SP)[4] to increase robustness to outliers and improve accuracy. We consider the problem of clustering faces given the number of categories for three face data sets, i.e., the Extended Yale Face Database B [7] (16128 images of 38 different subjects), the AR Face database [17] (4000 images of 126 different subjects) and the MUCT Face Database [19] (3755 images of 625 different subjects). All of them contain people under different illumination conditions. In addition, MUCT and AR include pose variations, and in the case of AR people use different outfits (see Fig.4 for some examples).
In our experiments we use the Parallel Spectral Clustering in Distributed Systems (PSCDS) [4] method as the base code for spectral clustering, but just employ a simple desktop machine. The different SFR methods are incorporated to PSCDS as a preprocessing stage as follows. First, each image is described by the Gist [21] holistic descriptor with scales of orientations and blocks. This produces a vector of dimensions. The use of Gist instead of the original images has consistently produced an improvement in accuracy in the range of [15%, 20%]. Secondly, all the descriptors of a dataset are combined forming an observation matrix , where is the total number of images in the specific dataset. The rank of is the number of expected clusters . Then, the SFR method under evaluation recovers a subspace of rank from . The matrix is then used in the pipeline of PSCDS to compute the distance matrix , considering five nearest neighbours per sample, followed by the spectral clustering.
We have considered LMAFIT, ROSL, ROSL+, FRADM and FRNys as representatives of the SFR approaches. Additionally, we also compared against stateoftheart clustering techniques such as the Robust Subspace Segmentation by LowRank Representation (LRR) [15] and the Smooth Representation Clustering (SMR) [10], specifically designed for clustering purposes. The results of our evaluation are presented in Table IV, including the average clustering error (ce); the base time, i.e., time taken by the specific SFR method; and the total time, i.e., base time plus the time taken by the PSCDS. For LRR and SMR the total time is that produced by the method.
When considering the YaleB and AR datasets FRADM obtains the lowest clustering errors, 2.7% and 6.65% respectively. Moreover, FRADM and FRNys present the best balance between accuracy and computational time for these datasets. MUCT, is the most challenging dataset with 625 classes, which is a very high rank in comparison to its matrix dimensions (). These conditions are beyond the recovery boundaries of SFR methods, and even though FRADM accuracy is comparable to that obtained by the top method, LRR. Furthermore, FRADM computational performance is more than 20 times faster than LRR for this case, supporting the good accuracyspeed tradeoff offered by the method.
PSCDS  LRR  SMR  LMAFIT  ROSL  ROSL+  FRADM  FRNys  
YaleB  ce (%)  18.7%  13.8%  28.4%  18.8%  20.1%  30.4%  2.7%  2.89% 
A=16128x5760  base time  5.17  64.8  351.6  2.4  274.6  7.6  6.8  0.58 
=38  total time  5.17  64.8  351.6  5.6  275.5  8.7  8.8  2.5 
AR  ce (%)  17.2%  36.8%  39.7%  6.70%  7.17%  46.0%  6.65%  7.17% 
A=4000x5760  base time  5.08  606.8  105.1  17.6  662.2  48.1  13.81  1.63 
=126  total time  5.08  606.8  105.1  21.05  665.1  49.9  16.91  3.83 
MUCT  ce (%)  55.3%  53.4%  55.8%  56.2%  56.3%  78.4%  55.7%  62.8% 
A=3755x5760  base time  76.7  3820  3995  101.2  17696  4890  85.5  67.2 
=625  total time  76.7  3820  3995  190.9  17771  4977  175.2  156.3 
In this paper we have proposed an efficient, stable and accurate technique, FRADM, to perform a robust decomposition of a corrupted matrix into its fixedrank and sparse components. To this end we have based our algorithm on a polar factorization on a product manifold , combining key tools from manifold optimization and fast projectors. We also proposed a fast SPD projector to speed up computation, along with a proof of its validity in this context. Additionally, Nystrom’s sampling techniques have been used to further accelerate the results, achieving a linear complexity. The resulting algorithm has been tested on synthetic cases and the challenging problems of robust photometric stereo and spectral clustering, proving to be as accurate and more efficient than stateoftheart approaches and paving the way towards largescale problems.
In this Appendix we shall show that the minimization subproblem (9), i.e.
(13) 
although highly nonconvex, converges geometrically to the global minimum when optimized via the proposed FixedRankOptFull method (Algorithm 4).
The FixedRankOptFull algorithm performs an alternating directions minimization (ADM) on each of the submanifolds , and (Algorithm 3). In each iteration it uses the algorithm FixedRankOptStep, described in Sec. 4.1, that performs a single step of the alternating directions minimization.
In Sec. 4.1 we provided the exact projectors on each of the submanifolds , and , and proved the validity of the ones corresponding to the Stiefel manifolds. For the case of the manifold, a careful analysis is required to prove its validity.
Given that , and considering and as the solutions of an OPP, then a unique solution in the SPD manifolds must exist. This solution is given in the following discussion, but we need some previous results.
Proof: Since , if , then and therefore , which is SPD since . A similar argument, but without any additional assumption on , since , shows that after minimizing with respect to .
Note that in Lemma 2 it is not neccesary that , only that it is invertible. Also, we conclude that .
, is in general not symmetric, although it can be written as a product of two SPD matrices, and therefore has positive eigenvalues. Even though from Lemma
2 we have that and are in , we cannot directly prove that , but we can do it passing to the limit inside Algorithm 4. When passing to the limit the sequences defined by , and , both conditions are simultaneously met, as in Lemma 3:Suppose that the FixedRankOptFull Algorithm converges to a fixed point , then , , and .
Proof: Since , and are solutions to their respective OPPs and have to be in their respective Stiefel manifolds. Then, by applying Lemma 2, both and are in SPD. Since , we have that:
(14)  
which is on SPD since it is a convex manifold. Then, by taking the square root of we have that .
Now, since the eigenvalues are continuous functions of the matrix entries, there exists such that all symmetric matrices in the open ball of radius centered at are contained in . Thus, if FixedRankOptFull converges, then there exists such that .
Let us now discuss the convergence of the FixedRankOptFull Algorithm. Given , then is the projector onto the column space of in . Note that , where . Then we have the following:
If , the FixedRankOptFull algorithm converges Qlinearly to a global minimum of (9) given by such that is the unique projection of onto . The convergence is Qlinear, in the sense that and .
Proof: For each , denote by the projectors as defined before. Then it is easy to proof, using the alternative definition of , that , where . Thus the sequence of subspaces is the same as that produced by the Orthogonal Iteration [8] for the computation of the first
eigenvalues and eigenvectors of the symmetric matrix
. The Orthogonal Iteration converges Qlinearly in the sense that , with the eigenvalues of . Since , we have that . By a similar argument .In our case, , with and is a perturbation matrix, then will be much smaller than and the error will be largely decreased in each iteration.
We would like to stress that although we do not provide an algebraic proof for due to its complexity, Lemma 3 along with the continuity of eigenvalues argument guarantee that when we are near an optimum. Starting with then for the first iteration we have , and according to Figure 1 very near the optimum, thus we can ensure that the whole sequence is in . This is not a complete proof, but Theorem 4 ensures global convergence despite the nature of . Thus at some point will be in , which is also shown to always occur in our extensive numerical experiments, even starting from random .
From few to many: illumination cone models for face recognition under variable lighting and pose.
Pattern Analysis and Machine Intelligence, IEEE Transactions on, 2001.
Comments
There are no comments yet.