Coordinate Descent And Proximal Algorithms For Precision Estimation

Why Optimization Structure Matters

The graphical lasso estimates a sparse precision matrix by balancing Gaussian likelihood fit against an off-diagonal $L_1$ penalty. Its objective is convex, but the positive-definite constraint and the nonsmooth penalty make the choice of optimization method consequential. A useful algorithm must do more than produce small changes in entries: it must respect the matrix domain, decrease the objective, and provide a meaningful check that it has approached a solution.

Two common approaches exploit different structures. Block coordinate descent works with a dual covariance-like matrix and turns each column update into a lasso problem. Proximal gradient works directly with the precision matrix: it takes a smooth likelihood step and soft-thresholds off-diagonal entries. The procedures are distinct, so their update equations and safeguards should not be mixed.

Definition: Let $S$ be a symmetric sample covariance matrix and let $\lambda\geq 0$. The graphical-lasso objective, using the convention that both symmetric off-diagonal entries are counted, is $$ F(\Theta)=-\log\det(\Theta)+\operatorname{tr}(S\Theta) +\lambda\sum_{i\ne j}|\Theta_{ij}|,\qquad \Theta\succ0. $$ Here $\Theta$ is symmetric positive definite. The diagonal entries are not penalized. The factor of two implicit in summing over $i\ne j$ is important when writing coordinate updates.

The likelihood part is smooth on the positive-definite cone. Writing $f(\Theta)=-\log\det(\Theta)+\operatorname{tr}(S\Theta)$, its gradient under the Frobenius inner product on symmetric matrices is $\nabla f(\Theta)=S-\Theta^{-1}$. The remaining term $g(\Theta)=\lambda\sum_{i\ne j}|\Theta_{ij}|$ is convex but not differentiable when an entry is zero.

Block Coordinate Descent Through the Dual

The dual formulation makes column-wise updates especially convenient. Introducing a symmetric matrix $W=S+U$, where $U_{ii}=0$ and $|U_{ij}|\leq\lambda$, gives the dual problem

\max_{W\succ0}\ \log\det(W) \quad\text{subject to}\quad W_{ii}=S_{ii},\qquad |W_{ij}-S_{ij}|\leq\lambda\quad(i\ne j).

An additive constant $p$ in the dual objective has been omitted. To see the relationship, express each penalty as $\lambda|\Theta_{ij}|=\max_{|U_{ij}|\leq\lambda}U_{ij}\Theta_{ij}$, with the symmetric entries treated consistently. For fixed feasible $W=S+U$, minimizing $\operatorname{tr}(W\Theta)-\log\det(\Theta)$ over $\Theta\succ0$ gives $\Theta=W^{-1}$: the derivative is $W-\Theta^{-1}$, and the minimized value is $p+\log\det(W)$. Thus the dual objective and constraints follow. At a solution, the primal precision estimate is the inverse of the dual matrix.

Now hold all entries of $W$ fixed except its $j$th row and column. Let $A=W_{-j,-j}$, let $w=W_{-j,j}$, and let $s=S_{-j,j}$. The diagonal entry remains $W_{jj}=S_{jj}$. The determinant identity gives $\det(W)=\det(A)(S_{jj}-w^\mathsf{T}A^{-1}w)$. Since $A$ is fixed, maximizing $\log\det(W)$ over this block is equivalent to minimizing $w^\mathsf{T}A^{-1}w$, subject to $|w-s|\leq\lambda$ componentwise and a positive Schur complement.

Result: If $A\succ0$ and the current dual matrix is feasible and positive definite, the $j$th block update can be computed by solving $$ \widehat\beta=\arg\min_{\beta\in\mathbb R^{p-1}} \left\{\frac12\beta^\mathsf{T}A\beta-s^\mathsf{T}\beta+\lambda\|\beta\|_1\right\}, \qquad w_{\mathrm{new}}=A\widehat\beta. $$ The updated block satisfies the dual box constraints and has a positive Schur complement.

For the derivation, set $\beta=A^{-1}w$, so that $w=A\beta$. The box constraint becomes $|A\beta-s|\leq\lambda$. The quadratic program $\min_\beta \frac12\beta^\mathsf{T}A\beta$ subject to that box constraint has Lagrange dual equal to the displayed lasso problem (the dual variable associated with $A\beta-s$ is $\beta$); equivalently, its optimality conditions are $A\beta-s+\lambda z=0$, where $z_k\in\partial|\beta_k|$. These conditions imply $|(A\beta-s)_k|\leq\lambda$, with equality and the appropriate sign when $\beta_k\ne0$. Finally, feasibility of the current positive-definite block supplies a feasible $w$ with $w^\mathsf{T}A^{-1}w<S_{jj}$. The minimizing $w_{\mathrm{new}}$ has no larger quadratic value, so its Schur complement is also positive. This proves that the block update remains in the positive-definite dual domain. $\square$

Worked Example: One Dual Column Update

Consider a two-variable update with $S_{jj}=4$, $A=(3)$, $s=(1.2)$, and $\lambda=0.5$. The one-dimensional lasso problem is $$ \min_{\beta}\ \frac12(3)\beta^2-1.2\beta+0.5|\beta|. $$

For a positive solution, its derivative is $3\beta-1.2+0.5$, giving $\widehat\beta=0.7/3=0.2333\ldots$. This is positive, so the assumed sign is consistent. The updated covariance-like entry is $w=A\widehat\beta=3(0.2333\ldots)=0.7$. It obeys the box constraint because $|w-s|=|0.7-1.2|=0.5=\lambda$.

The updated $2\times2$ dual block has diagonal entries $3$ and $4$, and off-diagonal entry $0.7$. Its determinant is $3(4)-0.7^2=11.51>0$, so it is positive definite. Its inverse is $$ \frac{1}{11.51} \begin{pmatrix}4&-0.7\\-0.7&3\end{pmatrix}. $$ This inverse is the corresponding primal precision estimate if this is the full matrix.

In practice, cycle through columns, solving each lasso subproblem with the currently updated $A$ and $s$, until the dual objective and the entries change sufficiently little. The current matrix must remain symmetric: a column update also updates its matching row. Because the subproblem is convex, a coordinate update optimizes that block with the others fixed; repeated cyclic updates target the dual optimum under the usual existence conditions for a positive-definite optimum.

Proximal Gradient in the Precision Matrix

A direct primal method separates the smooth likelihood $f$ from the penalty $g$. Given a current positive-definite iterate $\Theta_k$ and step size $\alpha_k>0$, first take a gradient step $Y_k=\Theta_k-\alpha_k(S-\Theta_k^{-1})$. Then apply the proximal map of $g$. On symmetric matrices with the Frobenius norm, the resulting entrywise update is

(\Theta_{k+1})_{ii}=(Y_k)_{ii},\qquad (\Theta_{k+1})_{ij}=\operatorname{sign}((Y_k)_{ij}) \max\{|(Y_k)_{ij}|-\alpha_k\lambda,0\},\quad i\ne j.

The threshold is $\alpha_k\lambda$, not $2\alpha_k\lambda$. Although the penalty counts both $(i,j)$ and $(j,i)$, the squared Frobenius distance also counts both entries. Minimizing over the shared symmetric value gives the stated threshold.

This formula is the unconstrained proximal step for the penalty, not a guarantee of positive definiteness. Soft-thresholding can move an eigenvalue through zero. A safeguarded implementation therefore accepts a candidate only if it is positive definite and satisfies the smooth majorization check

f(\Theta_{k+1})\leq f(\Theta_k) +\langle\nabla f(\Theta_k),\Theta_{k+1}-\Theta_k\rangle +\frac{1}{2\alpha_k}\|\Theta_{k+1}-\Theta_k\|_F^2.

If either check fails, reduce the trial step size (for example, multiply it by a fixed factor $0<\rho<1$) and recompute the candidate. The comparison uses $f$, not $F$, because the proximal minimization already accounts for $g$. The resulting accepted step satisfies $F(\Theta_{k+1})\leq F(\Theta_k)$, as shown below.

Worked Example: A Proximal Step and a Positive-Definiteness Check

Let $S=\begin{pmatrix}1&0.8\\0.8&1\end{pmatrix}$, $\Theta_k=I$, $\lambda=0.2$, and trial step $\alpha_k=0.5$. Since $\Theta_k^{-1}=I$, the gradient is $S-I=\begin{pmatrix}0&0.8\\0.8&0\end{pmatrix}$. The gradient step has diagonal entries $1$ and off-diagonal entries $-0.4$. Thresholding by $\alpha_k\lambda=0.1$ gives $\Theta_{k+1}=\begin{pmatrix}1&-0.3\\-0.3&1\end{pmatrix}$.

Its eigenvalues are $0.7$ and $1.3$, so the candidate is positive definite. Also $\|\Theta_{k+1}-\Theta_k\|_F^2=2(0.3)^2=0.18$, and the gradient inner product is $2(0.8)(-0.3)=-0.48$. The right side of the majorization check is $f(\Theta_k)-0.48+0.18=2-0.48+0.18=1.70$.

The candidate's smooth objective is $f(\Theta_{k+1})=-\log(0.91)+\operatorname{tr}(S\Theta_{k+1})$. The trace is $2-0.48=1.52$, so $f(\Theta_{k+1})\approx0.09431+1.52=1.61431\leq1.70$. Its penalty is $0.2(0.3+0.3)=0.12$, whereas the starting penalty is zero. Thus $F$ decreases from $2$ to approximately $1.73431$.

Descent and Convergence Checks

The descent guarantee follows directly from the proximal subproblem. Its minimizer obeys $g(\Theta_{k+1})+\langle\nabla f(\Theta_k),\Delta_k\rangle+\|\Delta_k\|_F^2/(2\alpha_k)\leq g(\Theta_k)$, where $\Delta_k=\Theta_{k+1}-\Theta_k$. Add the accepted majorization inequality for $f$ to both sides. The inner products cancel, yielding the stronger bound

F(\Theta_k)-F(\Theta_{k+1})\geq\frac{\|\Delta_k\|_F^2}{2\alpha_k}.

For a precise convergence statement, suppose the iterates remain in a compact subset of the positive-definite cone, and the backtracking procedure accepts a step satisfying the majorization check with $\alpha_k\in[\underline{\alpha},\overline{\alpha}]$, where $0<\underline{\alpha}\leq\overline{\alpha}<\infty$. These conditions can be enforced or verified in applications; compactness means eigenvalues are bounded above and bounded away from zero. Since $F$ is bounded below on that set, summing the descent bound gives

\sum_{k=0}^{\infty}\frac{\|\Delta_k\|_F^2}{2\alpha_k} \leq \frac{F(\Theta_0)-\inf_k F(\Theta_k)}{1}<\infty.

Using the upper bound $\alpha_k\leq\overline{\alpha}$, not the lower bound, now gives $\sum_k\|\Delta_k\|_F^2\leq2\overline{\alpha}\sum_k\|\Delta_k\|_F^2/(2\alpha_k)<\infty$. Hence $\|\Delta_k\|_F\to0$. This direction of the step-size bound is essential.

The proximal optimality condition is $0\in\nabla f(\Theta_k)+\Delta_k/\alpha_k+\partial g(\Theta_{k+1})$. Along any convergent subsequence of iterates, compactness gives a further subsequence for which $\alpha_k$ converges to some $\alpha_*\in[\underline{\alpha},\overline{\alpha}]$. Since $\Delta_k\to0$, continuity of $\nabla f$ on the positive-definite cone and the closedness of the convex subgradient relation imply that the subsequential limit $\Theta_*$ satisfies $0\in\nabla f(\Theta_*)+\partial g(\Theta_*)$. This is the first-order optimality condition for the convex objective, so every such cluster point is a minimizer. Because $-\log\det(\Theta)$ is strictly convex on the positive-definite cone and the other objective terms are convex, the minimizer is unique when it exists. If the iterates lie in the stated compact set, they therefore converge to that minimizer. $\square$

Worked Example: Checking a Candidate by Its KKT Residual

For $S=\begin{pmatrix}2&0\\0&1\end{pmatrix}$ and $\lambda=0.3$, consider $\Theta_*=\begin{pmatrix}0.5&0\\0&1\end{pmatrix}$. It is positive definite and $\Theta_*^{-1}=\begin{pmatrix}2&0\\0&1\end{pmatrix}=S$, so the smooth gradient is zero.

At zero off-diagonal entries, the subgradient condition for the symmetric off-diagonal coordinate is $|(S-\Theta_*^{-1})_{12}|\leq\lambda$. Its left side is zero, so the condition holds. The diagonal coordinates are unpenalized and their gradient entries are zero as well. Thus $0\in\nabla f(\Theta_*)+\partial g(\Theta_*)$, and convexity makes $\Theta_*$ a global minimizer.

A practical stopping check can combine a small relative objective decrease, a small relative change in $\Theta$, and a small stationarity residual. For proximal gradient, one convenient residual is $\|\Theta_k-\operatorname{prox}_{\alpha_k g}(\Theta_k-\alpha_k\nabla f(\Theta_k))\|_F/\alpha_k$. A small change alone is not enough: a very small step size can make the change tiny while the optimality residual remains substantial.

Numerical Safeguards and Common Pitfalls

Positive definiteness should be checked numerically, not inferred from an algorithm's label. For a symmetric candidate, an eigenvalue or Cholesky check can detect loss of positive definiteness; a small eigenvalue tolerance should reflect numerical precision and the scale of the matrix. In proximal gradient, reject a candidate if it is not safely positive definite, then reduce the step and recompute. Avoid silently replacing it with a non-symmetric matrix or continuing to evaluate $\log\det$ after the determinant has become nonpositive.

For dual block descent, maintain symmetry, fixed diagonal entries, and the box constraints. A lasso solver's convergence tolerance affects the accuracy of the block solution; checking the dual feasibility residual $\max_{i\ne j}(|W_{ij}-S_{ij}|-\lambda)_+$ is useful. Monitor the dual log-determinant as well as changes in $W$, and only invert the final positive-definite dual estimate to obtain $\Theta$.

Wrong: “The soft-thresholded matrix is symmetric, so it must be a valid precision matrix.” Correction: Symmetry does not imply positive definiteness; check eigenvalues or a Cholesky factorization.
Wrong: “The penalty counts both symmetric entries, so the threshold must be $2\alpha\lambda$.” Correction: The Frobenius distance also counts both entries; the resulting threshold is $\alpha\lambda$.
Wrong: “A decreasing objective and tiny entry changes prove convergence.” Correction: Also monitor a stationarity residual and use step-size and domain safeguards; tiny changes can arise from an excessively small step.
Wrong: “Summability of $\|\Delta_k\|_F^2/\alpha_k$ implies summability of $\|\Delta_k\|_F^2$ from a lower bound on $\alpha_k$.” Correction: The needed implication uses an upper bound on $\alpha_k$.

These algorithms address computation for a specified penalty value; the choice and calibration of that value are separate modeling decisions. The central computational principle is to pair each update with checks that match its mathematical domain: dual feasibility and positive definiteness for block updates, and positive definiteness plus majorization for proximal steps.

Check Your Understanding

Use these questions to check the update formulas and their safeguards.

Why does the dual block problem reduce to a lasso problem, and which matrix is inverted to recover the precision estimate?
Under the stated convention for the penalty, what is the off-diagonal soft-threshold in a proximal-gradient step? Explain the role of symmetry.
What two conditions should a backtracking implementation check before accepting a proximal candidate?
In the descent proof, which bound on the step sizes turns weighted square-summability into $\sum_k\|\Theta_{k+1}-\Theta_k\|_F^2<\infty$?
Why is a small successive-iterate difference not, by itself, a sufficient stopping criterion?