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:
- Can this code be made more efficient so that it can handle larger matrices (ideally with $n,m \approx 100$) much faster?
- 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)$)?
