Low Rank Matrix Completion

Matrix completion is the task of filling in the values of a partially observed matrix. The mathematical jargon is a bit heavy, but the code is relatively straightforward, and the result is well worth the effort!

Suppose we have a low-rank matrix $$X \in \mathbb{R}^{n \times m}$$ with missing values. Our goal is to find the complete matrix $$\tilde{X} \in \mathbb{R}^{n \times m}$$ such that $$\tilde{X}_{jk}=X_{jk}$$ for all non-missing entries $$(j,k)$$ in $$X$$.

Minimizing the Nuclear Norm

Rather than minimizing the rank of $$\tilde{X}$$ directly (which is NP-hard), we instead minimize its nuclear norm, $$|| \tilde{X} ||_* = \sum_j \sigma_j$$ (the sum over singular values). Using NMinimize works on a toy problem, but it is slow and scales poorly for larger matrices.

X=N@{{1,Indexed[x,{1,2}],3},{2,4,Indexed[x,{2,3}]},{Indexed[x,{3,1}],10,15}};
MatrixForm[X] -> MatrixForm[X/.Last@NMinimize[Total@SingularValueList@X,Variables@X]]
(* Actual result: {{1.,2.00028,3.},{2.,4.,5.9982},{4.99726,10.,15.}} *)
(* Ideal result: {{1.,2.,3.},{2.,4.,6.},{5.,10.,15.}} *)


Semidefinite Programming

A more efficient way to solve this problem is to recast it as semidefinite optimization as follows:

Minimize: trace($$U$$) + trace($$V$$)
Subject to: $$\tilde{X}_{jk}=X_{jk}$$ for non-missing $$(j,k)$$ in $$X$$
As well as: $$\begin{pmatrix} U & \tilde{X}\\ \tilde{X}^T & V \end{pmatrix} \succcurlyeq 0$$

The minimization is over all symmetric matrices $$U \in \mathbb{R}^{n \times n}$$ and $$V \in \mathbb{R}^{m \times m}$$ as well as all possible complete matrices $$\tilde{X}$$. The symbol "$$\succcurlyeq$$" denotes that the matrix is positive semidefinite. This method easily solves the toy problem above.

symmetricMatrix[var_,size_]:=Normal@SymmetrizedArray[{j_,k_}:>Indexed[var,{j,k}],{size,size},Symmetric[{1,2}]]
nuclearNormMinimization[X_]:=(
{n,m}=Dimensions@X;
U=symmetricMatrix[u,n];
V=symmetricMatrix[v,m];
mat=Transpose@Join[Transpose@Join[U,Transpose[X]],Transpose@Join[X,V]];
SemidefiniteOptimization[Tr[U]+Tr[V],VectorGreaterEqual[{mat,0},"SemidefiniteCone"],Variables@{X,U,V},Method->"SCS",Tolerance->10^-7]
)
X=N@{{1,Indexed[x,{1,2}],3},{2,4,Indexed[x,{2,3}]},{Indexed[x,{3,1}],10,15}};
MatrixForm[X] -> MatrixForm[X/.nuclearNormMinimization[X]]
(* Actual result: {{1.,2.,3.},{2.,4.,6.},{5.,10.,15.}} *)


Towards Better Efficiency and $$L_2$$ Regularization

The above code works for well when $$n,m \approx 10$$, but it slows down for larger matrices. The following 100x100 matrix has 10% of its values missing, and it requires a few seconds to complete.

lowRankMatrix[rank_,{n_,m_}]:=Sum[RandomReal[{-1,1},n]\[TensorProduct]RandomReal[{-1,1},m],{i,rank}];
SeedRandom[12345]
X=lowRankMatrix[5,{100,100}];
numMissing=Round[Times@@Dimensions@X/10];
pos=Transpose@{RandomInteger[{1,First@Dimensions@X},numMissing],
RandomInteger[{1,Last@Dimensions@X},numMissing]};
(X[[Sequence@@#]]=Indexed[x,#])&/@pos;
nuclearNormMinimization[X];//Timing
(* Timing: About 3 seconds *)


In addition, $$\tilde{X}$$ must exactly equal $$X$$ at the measured values, but ideally this constraint would be relaxed to account for potentially noisy data. Thus, I am interested in the following two modifications:

1. Can this code be made more efficient so that it can handle larger matrices (ideally with $$n,m \approx 100$$) much faster?
2. Can the minimization be modified to include an $$L_2$$ penalty (i.e. minimizing $$\text{trace(U)}+\text{trace(V)}+\sum_{jk} (\tilde{X}_{jk}-X_{jk})^2$$, where the sum is over all non-missing entries $$(j,k)$$)?

Robust PCA seems to handle this problem well, but the value for missing data might be tricky.

For Robust PCA the optimization problems is $$\min {|| L ||}_{*} + \lambda {|| S ||}_{1}$$ s.t. $$M = L + S$$, where $$M$$ is an input matrix, $$L$$ is a low rank matrix and $$S$$ is a sparse matrix. $${||~||}_{*}$$ stands for nuclear norm and $${||~||}_{1}$$ is $$L_1$$ norm.

lowRankMatrix[rank_,{n_,m_}]:=Sum[RandomReal[{-1,1},n]\[TensorProduct]RandomReal[{-1,1},m],{i,rank}];
SeedRandom[12345]
Y = X = lowRankMatrix[5,{40,60}];
numMissing=Round[Times@@Dimensions@X/10];
pos=Transpose@{RandomInteger[{1,First@Dimensions@X},numMissing],
RandomInteger[{1,Last@Dimensions@X},numMissing]};
X === Y
(X[[Sequence@@#]]=Indexed[x,#])&/@pos;
X = X /. Thread[Variables[X] -> -10.0] ;
X === Y
MATRIX = X ;
MU = 1/4*1/Norm[MATRIX, 1]*Apply[Times, Dimensions[MATRIX]] ;
LAMBDA = 1/Sqrt[N[Max[Dimensions[MATRIX]]]] ;
FACTOR = 10.0^-12 ;
TOLERANCE = FACTOR*Norm[MATRIX, "Frobenius"] ;
LIMIT = 10000 ;
DATA = RPCA[MATRIX, MU, LAMBDA, TOLERANCE, LIMIT] ;
DATA//Last
MatrixRank[First[DATA]]
MinMax[Abs@Chop[First[DATA]-Y]]


Note, this implementation doesn't have adaptive rank truncation, hyperparameters might not be optimal for your problem. Code for RPCA:

(* Rachel Thomas, Computational Linear Algebra for Coders, https://github.com/fastai/numerical-linear-algebra, nb2 & nb3 *)
ClearAll[RandomizedRange] ;
RandomizedRange[matrix_?MatrixQ, size_Integer, iteration_Integer, seed_:1] := Block[
{projection, transpose},
SeedRandom[seed] ;
projection = RandomVariate[NormalDistribution[], {Last[Dimensions[matrix]], size}] ;
transpose = Transpose[matrix] ;
Do[
projection = Transpose[First[QRDecomposition[Dot[matrix, projection]]]] ;
projection = Transpose[First[QRDecomposition[Dot[transpose, projection]]]],
iteration
] ;
Transpose[First[QRDecomposition[Dot[matrix, projection]]]]
] ;
ClearAll[RandomizedSingularValueDecomposition] ;
RandomizedSingularValueDecomposition::usage = "
RandomizedSingularValueDecomposition[m, k] -- gives the randomized singular value decomposition associated with the k largest singular values of m
" ;
Options[RandomizedSingularValueDecomposition] = {
"BufferSize" -> 10,
"NumberOfIterations" -> 10,
"RandomSeed" -> 1
} ;
RandomizedSingularValueDecomposition[matrix_?MatrixQ, size_Integer, options:OptionsPattern[]] := Block[
{buffer, iteration, seed, projection, u, s, v},
buffer = OptionValue["BufferSize"] ;
iteration = OptionValue["NumberOfIterations"] ;
seed = OptionValue["RandomSeed"] ;
projection = RandomizedRange[matrix, size + buffer, iteration, seed] ;
{u, s, v} = SingularValueDecomposition[Dot[Transpose[projection], matrix], size] ;
u = Dot[projection, u] ;
{u, s, v}
] ;
(* https://github.com/dganguli/robust-pca *)
ClearAll[RPCAShrink] ;
RPCAShrink[matrix_, tau_] := Sign[matrix]*(Ramp[Abs[matrix] - tau]) ;
ClearAll[RPCAThreshold] ;
RPCAThreshold[matrix_, tau_, svd_, rnk_] := Block[
{u, s, v},
{u, s, v} = svd[matrix, rnk] ;
Dot[u, Dot[RPCAShrink[s, tau], Transpose[v]]]
] ;
ClearAll[RPCA] ;
RPCA::usage = "
RPCA[matrix, mu, lambda, tolerance, limit] -- perform RPCA decomposition of matrix, returns {low rank, sparse, {number of iterations, error}}
" ;
Options[RPCA] = {
"RNK" -> Full,
"SVD" -> SingularValueDecomposition
} ;
(* hyper-parameters *)
(* mu = 1/4*1/Norm[matrix, 1]*Apply[Times, Dimensions[matrix]] ; *)
(* lambda = 1/Sqrt[N[Max[Dimensions[matrix]]]] ; *)
(* tolerance = 10.0^-8*Norm[matrix, "Frobenius"] ; *)
RPCA[matrix_?MatrixQ, mu_Real, lambda_Real, tolerance_Real, limit_Integer, options:OptionsPattern[]] := Catch[
Block[
{m, n, rnk, svd, inverse, count, error, sk, yk, lk},
{m, n} = Dimensions[matrix] ;
rnk = OptionValue["RNK"] ;
rnk = Which[
SameQ[rnk, Full],
Min[{m, n}],
MatchQ[rnk, rnk_Integer?Positive /; rnk <= Min[{m, n}]],
rnk,
Throw[$$Failed], True, Throw[$$Failed]
] ;
svd = OptionValue["SVD"] ;
If[Not[MemberQ[{SingularValueDecomposition, RandomizedSingularValueDecomposition}, svd]], Throw[\$Failed]] ;
inverse = 1.0/mu ;
count = 0 ;
sk = yk = lk = ConstantArray[0.0, Dimensions[matrix]] ;
While[
count < limit,
lk = RPCAThreshold[matrix - sk + inverse*yk, inverse, svd, rnk] ;
sk = RPCAShrink[matrix - lk + inverse*yk, inverse*lambda] ;
error = matrix - lk - sk ;
yk = yk + mu*error ;
error = Norm[error, "Frobenius"] ;
count++ ;
If[error < tolerance, Break[]] ;
] ;
{lk, sk, {count, error}}
]
] ;
RPCA[matrix_?MatrixQ, factor_Real, limit_Integer, options:OptionsPattern[]] := RPCA[
matrix,
1/4*1/Norm[matrix, 1]*Apply[Times, Dimensions[matrix]],
1/Sqrt[N[Max[Dimensions[matrix]]]],
factor*Norm[matrix, "Frobenius"],
limit,
options
] ;

• Seems like a powerful approach, but how would you robustly fill in missing values? For example, after the update step for lk, could you replace all its values with those of matrix at the non-missing positions? In other words, you would force the L matrix to exactly match the input matrix at all the measured values. Would this disrupt the algorithm? (After playing around with this, it seems to work well, with lambda=0 being the optimal parameter for my system.)
– Tal
Jun 18, 2021 at 3:20
• @Tal, here replacing missing values with -10 worked, but it is hard to say how to deal with the general case. It seems RPCA works best when missing values are viewed as grain noise, i.e. in some sense are far from true values. Perhaps a different matrix decomposition can perform better in your case, see this or this
– I.M.
Jun 19, 2021 at 3:10

Here is an approach I recently learned about that explicitly takes into account the missing values during the update stage. Provided all the entire of $$X$$ are positive, it finds a complete matrix $$\tilde{X}$$ that can be factorized as $$\tilde{X} = W H$$ where $$\tilde{W} \in \mathbb{R}^{n \times k}$$, $$\tilde{H} \in \mathbb{R}^{k \times m}$$, and $$k$$ is the rank of the factorization.

The update rules are very simple. Initialize $$H$$ and $$W$$ to be positive matrices. In each iteration, we perform entry-wise multiplication and division:

$$H \to H*\frac{W^T X}{W^T (W H)}$$
$$W \to W*\frac{X H^T}{(W H) H^T}$$

For all positions in $$X$$ with missing values, replace those positions in the product $$(W H)$$ within the denominators of the update rules.

matrixFactorize[Xraw_?MatrixQ, rank_] := Block[{posMissing, min, X, W, H},
posMissing = Position[Xraw, _Missing];
(* Ensures that all matrix elements are positive *)
min = Min@DeleteMissing@Flatten@Xraw;
X = Xraw - min;
X = ReplacePart[X, posMissing -> 0];
(* Initialize W and H *)
SeedRandom[12345];
W = RandomReal[{0, 1}, {Length@X, rank}];
H = RandomReal[{0, 1}, {rank, Length@First@X}];
current = RootMeanSquare@Flatten@Delete[W.H - X, posMissing];
(* Update the W and H matrices *)
Do[
H = H*Transpose[W].X/(Transpose[W].ReplacePart[W.H, posMissing -> 0]);
W = W*X.Transpose[H]/(ReplacePart[W.H, posMissing -> 0].Transpose[H]);
,{iter,100}];
W.H + min
]


This algorithm is very fast, with the example below filling in a 100x100 matrix with 10% missing values. However, from general experimentation, it appears less accurate at filling in missing values than nuclear norm minimization.

SeedRandom[12345]
Xcomplete = X = lowRankMatrix[5, {100, 100}];
numMissing = Round[Times @@ Dimensions@X/10];
X = ReplacePart[X, Transpose@{RandomInteger[{1, First@Dimensions@X}, numMissing], RandomInteger[{1, Last@Dimensions@X}, numMissing]} -> Missing[]];
RootMeanSquare@Flatten[Xcomplete - matrixFactorize[X, 5]] // Timing
(* Timing: 0.3 seconds *)
(* Root Mean Square Error: 0.3 *)