2.6 Low-rank matrix completion

Motivation

A partially observed table contains two different kinds of information: measured values and locations where no value has been measured. Matrix completion uses relationships across rows and columns to estimate those missing values. A low-rank model expresses these relationships through a small number of shared latent factors. The same factors connect observed and unobserved entries, allowing information from one part of the table to influence another.

The missingness pattern is as important as the number of measurements. Three entries can determine a small rank-one table, while another observation pattern can leave many possible completions. We will first specify exactly which observations enter the loss, then use a convex singular-value penalty to encourage a simpler fitted matrix. The resulting algorithm combines a gradient step on observed residuals with singular-value thresholding across the whole matrix.

Definitions

Let $P\in\mathbb R^{n\times d}$ store the available data, and let $\Omega\subseteq\{1,\ldots,n\}\times\{1,\ldots,d\}$ record observed positions. Entries of $P$ outside $\Omega$ are storage choices, not targets for the fit. An observed zero remains a measurement and belongs to $\Omega$.

Definition

Observation projection and matrix norms

The observation operator keeps precisely the observed coordinates:

$$(\mathcal P_\Omega X)_{ij}= \begin{cases}X_{ij},&(i,j)\in\Omega,\\0,&(i,j)\notin\Omega.\end{cases} \qquad \mathcal P_{\Omega^c}=I-\mathcal P_\Omega.$$

For matrices of the same shape, the Frobenius inner product and norm are $\langle X,Y\rangle_F=\sum_{ij}X_{ij}Y_{ij}$ and $\|X\|_F^2=\sum_{ij}X_{ij}^2$. If the singular values of $X$ are $\sigma_1,\ldots,\sigma_p$, where $p=\min(n,d)$, its nuclear norm is

$$\|X\|_*=\sum_{j=1}^p\sigma_j,\qquad \|X\|_F^2=\sum_{j=1}^p\sigma_j^2.$$

The two norms summarize different aspects of the singular spectrum. The Frobenius norm measures total squared size, while the nuclear norm adds singular magnitudes. The rank instead counts positive singular values. Penalizing the nuclear norm can set singular values to zero through a thresholding step, without explicitly searching over discrete rank values.

Observation geometry and low-rank factors

Theorem

The observation mask is an orthogonal projection

With the Frobenius inner product, $\mathcal P_\Omega$ is self-adjoint and idempotent. In particular,

$$\|X\|_F^2=\|\mathcal P_\Omega X\|_F^2+ \|\mathcal P_{\Omega^c}X\|_F^2,\qquad \|\mathcal P_\Omega X\|_F\le\|X\|_F.$$

Proof

Projection onto the measured coordinates

Both $\langle\mathcal P_\Omega X,Y\rangle_F$ and $\langle X,\mathcal P_\Omega Y\rangle_F$ equal $\sum_{(i,j)\in\Omega}X_{ij}Y_{ij}$, proving self-adjointness. Applying the same mask twice changes nothing, proving idempotence. The observed and unobserved parts have disjoint support, so their inner product is zero. Pythagoras gives the norm identity and its inequality.

A rank-at-most-$r$ model can be written $X=LR^T$, where $L\in\mathbb R^{n\times r}$ and $R\in\mathbb R^{d\times r}$. Entry $X_{ij}$ is the inner product of row $i$ of $L$ and row $j$ of $R$. These shared coordinates make a missing entry depend on the same factors that explain measured entries. Searching directly for the factors produces a nonconvex problem because they multiply each other.

The factors themselves are not uniquely named coordinates. For any invertible $Q\in\mathbb R^{r\times r}$, replacing them by $LQ$ and $RQ^{-T}$ leaves $LR^T$ unchanged. Even the completed matrix may be undetermined by the observations. To discuss recovery, distinguish this coordinate freedom from whether the missing entries have a unique value. Low rank supplies a structural assumption; the observation pattern must still contain enough connections to determine the desired entries.

Worked example

1. A determined entry and an ambiguous observation pattern

Suppose the observed entries give

$$X=\begin{bmatrix}2&b\\3&6\end{bmatrix}.$$

A rank-one completion requires $\det X=12-3b=0$, so $b=4$ uniquely. The completed matrix factors as $(2,3)^T(1,2)$.

If only the two diagonal entries are observed, all matrices

$$X_t=\begin{bmatrix}2&t\\12/t&6\end{bmatrix},\qquad t\ne0,$$

have rank one and match the observations. For instance, $t=2$ and $t=4$ give different missing values. A low-rank assumption alone cannot select which is the underlying table. A regularizer can express a preference among fits, but that preference does not add an observation.

A loss on observed entries

For measurements that may contain noise, consider the objective

$$\Phi(X)=f(X)+\lambda\|X\|_*,\qquad f(X)=\frac12\|\mathcal P_\Omega(X-P)\|_F^2,\qquad\lambda>0.$$

Only measured residuals enter $f$. Replacing the unused entries of $P$ by other numbers does not change the objective. The nuclear penalty acts on the whole candidate matrix and connects its coordinates through singular directions. Both terms are convex, although the full objective need not be strictly convex for an incomplete observation pattern.

The parameter controls a fit-and-complexity compromise. Larger values place more weight on shrinking singular magnitudes, often removing more singular directions. The fitted observed entries can differ from the measurements. For exact measurements, the separate constrained model $\min\|X\|_*$ subject to $\mathcal P_\Omega X=\mathcal P_\Omega P$ expresses exact agreement. Here we develop the displayed penalized objective, whose observed residuals remain part of the optimization.

The nuclear norm is a convex matrix norm. One useful justification is its dual description as the largest inner product with a matrix of spectral norm at most one. The singular-value trace inequality bounds that inner product by the sum of singular values, and aligned singular directions attain the bound. A maximum of linear functions is convex. This connects the convexity of the penalty to the same singular-vector geometry used below.

Singular-value thresholding

Definition

Soft thresholding of singular values

For $Y=U\operatorname{diag}(s_j)V^T$ and $\tau\ge0$, define

$$\operatorname{SVT}_\tau(Y) =U\operatorname{diag}\big((s_j-\tau)_+\big)V^T, \qquad (a)_+=\max(a,0).$$

Theorem

The nuclear-norm proximal problem

For any fully specified matrix $Y$ and $\tau\ge0$, the unique minimizer of

$$\frac12\|X-Y\|_F^2+\tau\|X\|_*$$

is $\operatorname{SVT}_\tau(Y)$.

Proof

Reducing the matrix problem to scalar thresholds

Let the ordered singular values of $X$ and $Y$ be $t_j$ and $s_j$. The von Neumann trace inequality states $\langle X,Y\rangle_F\le\sum_jt_js_j$. Expanding the squared distance therefore bounds the objective below by

$$\sum_j\left[\frac12(t_j-s_j)^2+\tau t_j\right],\qquad t_j\ge0.$$

Each summand is minimized at $t_j=(s_j-\tau)_+$. This operation preserves their order. Choosing the singular vectors of $X$ to align with those of $Y$ makes the trace inequality an equality and attains every scalar minimum. The squared-distance term is strictly convex, so the minimizing matrix is unique, even when a choice of singular-vector basis is not unique.

This threshold subtracts the same amount from each positive singular value and sets singular values at or below the threshold to zero. A rank truncation instead keeps the largest selected singular values unchanged. Both can produce lower-rank matrices, but they solve different objectives. At a repeated singular value, a numerical SVD can return a different orthonormal basis for its subspace; applying the same scalar threshold throughout that subspace gives the same reconstructed matrix.

Worked example

2. A complete-observation threshold calculation

Take $Y=\operatorname{diag}(3,1)$ and $\tau=1/2$. Its singular directions are the coordinate axes, so the solution is $X=\operatorname{diag}(5/2,1/2)$. Its squared-distance contribution is $\tfrac12(1/4+1/4)=1/4$, and its penalty is $\tfrac12(5/2+1/2)=3/2$. The total objective is $7/4$. Keeping $Y$ unchanged gives objective $2$, so exact data agreement is not the minimizer of this penalized problem. The threshold has shrunk both observed singular magnitudes.

A masked proximal-gradient step

The gradient of the observed loss is

$$\nabla f(X)=\mathcal P_\Omega(X-P),\qquad \|\nabla f(X)-\nabla f(Y)\|_F\le\|X-Y\|_F.$$

The inequality follows from the projection theorem, so one is a Lipschitz constant for the gradient. A step size $0<\eta\le1$ gives the update

$$X^+=\operatorname{SVT}_{\eta\lambda} \left(X-\eta\mathcal P_\Omega(X-P)\right).$$

For unit step size this becomes

$$X^+=\operatorname{SVT}_\lambda \left(\mathcal P_\Omega P+\mathcal P_{\Omega^c}X\right).$$

The forward point uses measurements at observed positions and the current estimate elsewhere. Thresholding then acts on that entire matrix. The output can change observed entries as well as fill missing ones, consistently with the objective's balance between observed residual and nuclear penalty. The unobserved coordinates survive in the next forward point and influence the next SVD; they are not repeatedly treated as measured zeros.

To justify the step, expand the quadratic loss at $X$ with increment $D$:

$$f(X+D)=f(X)+\langle\nabla f(X),D\rangle_F +\tfrac12\|\mathcal P_\Omega D\|_F^2 \le f(X)+\langle\nabla f(X),D\rangle_F+\tfrac12\|D\|_F^2.$$

Adding the nuclear penalty creates a quadratic upper model equal to $\Phi$ at $X$. Its minimizer is the unit-step update, by the proximal theorem. Minimizing this upper model cannot increase $\Phi$. For a smaller step, replacing its quadratic coefficient by $1/(2\eta)$ gives another valid upper model. The step-size restriction therefore comes directly from the mask's norm, rather than from an arbitrary tuning rule.

Penalty size and an exact zero-solution test

There is a useful scale for the regularization parameter. Put $Y=\mathcal P_\Omega P$ and write $\|Y\|_2$ for its largest singular value. The trace inequality gives

$$\Phi(X)-\Phi(0) =\tfrac12\|\mathcal P_\Omega X\|_F^2-\langle Y,X\rangle_F +\lambda\|X\|_* \ge\tfrac12\|\mathcal P_\Omega X\|_F^2 +(\lambda-\|Y\|_2)\|X\|_*.$$

Consequently, if the parameter is at least this largest singular value, zero is a global minimizer. This gives an exact explanation for an iteration that removes every singular direction. The parameter has assigned enough cost to matrix magnitude that any possible improvement in observed fit is outweighed. The threshold scale is computed from the observed table and its mask, not from an imagined fully observed matrix. It also supplies a useful endpoint when examining a sequence of parameter choices, with smaller parameters allowing increasingly substantial fitted structure.

What a fixed point certifies

The gradient alone does not describe optimality for a nonsmooth penalized objective. A subgradient $G$ of the nuclear norm at $X$ means that $\|Z\|_*\ge\|X\|_*+\langle G,Z-X\rangle_F$ for every $Z$. The unit-step proximal update satisfies

$$\frac{X-X^+-\nabla f(X)}{\lambda}\in\partial\|X^+\|_*.$$

If an update returns the same matrix, this supporting inequality balances the observed gradient against the penalty. Adding it to the ordinary first-order inequality for the convex function $f$ proves that every candidate matrix has objective at least the current one. A fixed point is therefore a global minimizer of the stated convex objective. Observed residuals may remain nonzero at this point because their gradient is balanced by a nuclear-norm subgradient.

This argument motivates checking the size of the proximal update as well as the objective. A small change tests approximate stationarity of the complete penalized problem. It does not require the observed residual gradient to vanish, nor does it establish how accurate unobserved predictions are.

Worked example

3. One step fills an unobserved entry

Use the three observed entries of the first example, store $P=\begin{bmatrix}2&0\\3&6\end{bmatrix}$, initialize $X_0=0$, and take $\lambda=1$. The missing upper-right position is excluded from $\Omega$. The forward point is $Y=P$, whose singular values are $(\sqrt{73}+5)/2$ and $(\sqrt{73}-5)/2$. Both exceed one. Its orthogonal polar factor is

$$UV^T=\frac1{\sqrt{73}}\begin{bmatrix}8&-3\\3&8\end{bmatrix},\qquad X_1=\begin{bmatrix} 2-8/\sqrt{73}&3/\sqrt{73}\\ 3-3/\sqrt{73}&6-8/\sqrt{73} \end{bmatrix}.$$

Thus the missing estimate becomes about $0.351123$, even though it was zero in the initial storage table. The observed estimates also move. The objective falls from $49/2=24.5$ to $\sqrt{73}-2+137/146\approx7.482360$. The exact rank-one completion has missing value $4$. The first penalized iteration gives an intermediate estimate under its fit-and-penalty criterion.

Executable calculation

The Boolean mask records observation status independently of numeric values. The following implementation uses unit steps and the objective just derived.

import numpy as np

observed = np.array([[True, False], [True, True]])
P = np.array([[2., 0.], [3., 6.]])
lam = 1.

def svt(Y, threshold):
    U, s, Vt = np.linalg.svd(Y, full_matrices=False)
    return (U * np.maximum(s - threshold, 0.)) @ Vt

def objective(X):
    residual = np.where(observed, X - P, 0.)
    return 0.5 * np.sum(residual**2) + lam * np.linalg.svd(
        X, compute_uv=False
    ).sum()

X = np.zeros_like(P)
first = svt(P, lam)
history = [objective(X)]
for _ in range(100):
    gradient = np.where(observed, X - P, 0.)
    X = svt(X - gradient, lam)
    history.append(objective(X))

assert np.all(np.diff(history) <= 1e-10)
print(np.round(first, 6))
print(round(history[0], 6), round(history[-1], 6))

The first matrix printed is approximately $\begin{bmatrix}1.063671&0.351123\\2.648877&5.063671\end{bmatrix}$. The objective values printed are $24.5$ and approximately $7.073338$. The finite iteration count gives a reproducible numerical example. For a larger application, report the objective and a measure of iterate change when deciding when to stop. A low fitting loss alone cannot measure accuracy at positions whose values were never supplied.

Scaling the parameter and assessing the fit

The numerical size of the parameter depends on how values and the loss are scaled. If all measurements are multiplied by a positive constant $c$, scale the parameter by the same constant. For a candidate scaled in the same way, both the residual term and the penalty are then multiplied by $c^2$, so the ordering of all candidate objectives is unchanged. This identifies how an equivalent model behaves after a change of measurement units, rather than assigning an absolute interpretation to one parameter value.

Similarly, replacing the summed observed squared loss by its average requires dividing the penalty coefficient by the number of observations to obtain the same optimizer. Keeping the coefficient unchanged while averaging the loss instead increases the penalty's relative influence. State the loss normalization before comparing parameter values across implementations or tables with different amounts of observed data.

To assess predictions, reserve some genuinely measured entries for validation and construct the fitting mask from the remaining entries. Evaluate predicted values on the held-out positions without including those positions in gradient steps. This measures predictions of observations not used for fitting. It also changes the fitting observation pattern, so the split should retain the row and column information needed by the intended model. Optimization diagnostics and prediction diagnostics answer distinct questions: one checks whether the selected objective has been minimized, while the other checks the usefulness of its resulting estimates.

Common mistakes

Common mistake

Using zero values to infer the observation mask

An observed zero contributes its squared residual just like any other measurement. A missing entry contributes no residual, regardless of the number stored there. Construct the mask from observation status. Also keep the order of the derived update: make the gradient forward point, then apply singular-value thresholding. This preserves the stated penalized objective and its step-size argument.

Summary

The observation operator is an orthogonal projection, which determines both the loss and its gradient. Low-rank structure couples missing entries, while identifiability depends on where measurements occur. The nuclear norm gives a convex singular-value penalty, and its proximal map soft-thresholds the singular values. Combining that map with the observed gradient yields an algorithm whose unit step has a direct quadratic-upper-bound justification.

Exercises

Checkpoint

1. A rank-at-most-one table has observed entries $X_{11}=0$, $X_{21}=2$, and $X_{22}=4$. Determine the missing $X_{12}$. Explain why the observed zero belongs in the mask.

Use the determinant and the observation definition.

Checkpoint

2. For full observations $Y=\operatorname{diag}(5,2)$ and penalty parameter $\tau=3$, find the nuclear proximal solution, its rank, and the minimized objective.

Compute both the squared residual term and the nuclear penalty.

Checkpoint

3. In the masked example, change the penalty to $\lambda=7$ and start at zero. What is the first proximal-gradient update? Does the objective impose exact agreement with the three observed entries?

Compare the threshold with the larger displayed singular value.

Solutions

Solution · 1. A zero measurement supplies information

The determinant is $0\cdot4-2X_{12}$, so rank at most one forces $X_{12}=0$. The second row is nonzero, so the completed table has rank one. Omitting $X_{11}=0$ would remove an actual equation and permit other rank-one completions. Its zero value does not make it missing.

Solution · 2. Thresholded singular magnitudes

The solution is $\operatorname{diag}(2,0)$, of rank one. Its squared-distance term is $\tfrac12(3^2+2^2)=13/2$, and the penalty is $3\cdot2=6$. The minimized objective is $25/2=12.5$.

Solution · 3. A penalty that removes every initial singular value

The largest forward-point singular value is $(\sqrt{73}+5)/2<7$. Thresholding at seven gives the zero matrix for the first update. The model uses an observed-residual penalty, so it permits disagreement with measured entries; it does not impose their equality as a constraint.

Related notes

• SVD and principal components
• Least squares and residuals
• Gradient descent on quadratics
• Course reading