Setup

Recall the $L^2$ regularized least squares problem (i.e., ridge regression),

\begin{align} u_{\star} &:= \text{argmin}_{u \in \mathbb{R}^d} J(u) \tag{1} \newline J(u) &:= \frac{1}{2} \lVert y - Gu\rVert^2_{\Sigma} + \frac{1}{2}\lVert u-m \rVert^2_{C} \end{align}

where $y \in \R^n$ and $m \in \R^d$ are fixed vectors, $G \in \R^{n \times d}$ is the linear forward model, and $\Sigma \in \R^{n \times n}$ and $C \in \R^{d \times d}$ are positive definite matrices. In (1) we have used the following notation for inner products and norms weighted by a positive definite matrix $C$: \begin{align} \langle v, v^\prime \rangle_C &:= \langle C^{-1}v, v^\prime\rangle = v^\top C^{-1}v^\prime \newline \lVert v \rVert_C^2 &:= \langle v, v\rangle_C. \end{align}

A solution to (1) can be viewed as a maximum a posteriori (MAP) estimate under the Bayesian model \begin{align} y|u &\sim \mathcal{N}(Gu, \Sigma) \tag{3} \newline u &\sim \mathcal{N}(m, C). \tag{4} \end{align}

I discuss this problem in depth in my post, on the linear Gaussian model, where I show that the solution can be written as the following equivalent expressions, \begin{align} u_{\star} &= \left(G^\top \Sigma^{-1}G + C^{-1}\right)^{-1}\left(G^\top \Sigma^{-1}y + C^{-1}m \right) \tag{5} \newline u_{\star} &= m + CG^\top \left(GCG^\top + \Sigma \right)^{-1}(y-Gm). \tag{6} \end{align}

In this post we will consider a generalization of problem (1), where $C$ is only required to be positive semidefinite (not strictly positive definite). In particular, this means that the inverse $C^{-1}$ may not exist. It is interesting to note that, although (5) depends on $C^{-1}$, the inverse does not appear anywhere in the equivalent expression (6). This leads to the natural question of whether (6) still provides a valid solution to the optimization problem even when $C$ is singular. The first step in investigating this question will be providing a suitable generalization of the objective function in (1), since the expression $\lVert u - m\rVert^2_{C} = \lVert C^{-1/2}(u - m)\rVert^2$ is not well-defined when $C$ is singular.

Generalizing the Optimization Problem

Our goal is now to extend (1) to the case where $C$ need not be strictly positive definite. We start by generalizing (1) through the definition of a constrained optimization problem on an extended parameter space. We then provide two approaches to formulate this problem in an unconstrained fashion, each of which leads to equivalent closed-from solutions. We will find it convenient to define notation for the “centered” quantities \begin{align} &u^\prime := u - m, &&y^\prime := y - Gm \tag{7} \end{align} so that $y - Gu = (y - Gm) - G(u-m) = y^\prime - Gu^\prime. \tag{8}$

Constrained Formulation

We start by considering the constrained formulation. This section follows the approach in section 8.1.2 of (Sanz-Alonso et al., 2023). To start, note that when $C$ is positive definite, we have $\lVert u - m\rVert^2_C = \lVert u^\prime \rVert^2_C = \langle C^{-1}u^\prime, u^\prime\rangle = \langle b, u^\prime \rangle, \tag{9}$ where $b \in \mathbb{R}^d$ solves $Cb = u^\prime. \tag{10}$ When $C$ is invertible, the unique $b$ solving (10) is simply found by multiplying both sides by $C^{-1}$. When $C$ is not invertible, then there may be zero or infinitely many such solutions. As long as there is at least one solution to (10) we can give meaning to expression (9) by picking a particular $b$ that solves (10). Now that we have introduced a new variable $b$, we might consider generalizing (1) by jointly optimizing over $(u,b)$ subject to the linear constraint $Cb = u^\prime$.

Constrained Joint Optimization.
\begin{align} (u_{\star}, b_{\star}) &\in \text{argmin}_{(u,b) \in \mathcal{S}} J(u,b) \tag{11} \newline J(u,b) &:= \frac{1}{2} \lVert y - Gu\rVert^2_{\Sigma} + \frac{1}{2}\langle b, u-m \rangle \newline \mathcal{S} &:= \left\{(u,b) : Cb = u-m \right\}. \end{align}

Note that if $C$ is positive definite, then (11) reduces to (1). If the problem is solved for $(u_{\star}, b_{\star})$, then the desired solution can be obtained by extracting $u_{\star}$ and discarding the nuisance parameter $b_{\star}$.

Lagrange Multipliers

Observe that (11) can be viewed as optimizing $J(u,b)$ subject to the constraint $g(u,b) = 0$, where $g(u,b) = Cb - (u-m)$. This is a typical setup for Lagrange multipliers. We therefore introduce the Lagrange multiplier $\lambda \in \mathbb{R}^d$, which allows us to cast the constrained optimization over $(u,b)$ as an unconstrained optimization over $(u,b,\lambda)$.

Lagrange Multiplier Formulation.
\begin{align} (u_{\star}, b_{\star}, \lambda_{\star}) &\in \text{argmin}_{u,b,\lambda} J(u,b,\lambda) \tag{12} \newline J(u,b,\lambda) &:= \frac{1}{2} \lVert y - Gu\rVert^2_{\Sigma} + \frac{1}{2}\langle b, u-m \rangle + \langle \lambda, Cb - u + m \rangle. \end{align}

We have succeeding in converting the problem to one that can be solved by analytical means; namely, by solving the system of equations $\nabla_{u} J = \nabla_{b} J = \nabla_{\lambda} J = 0. \tag{13}$

This derivation is provided in the appendix, and the result is summarized below.

Solution of Lagrange Multiplier Formulation.
A solution $(u_{\star}, b_{\star}, \lambda_{\star})$ of (12), projected onto the $u$-component, is given by $u_{\star} = m + Cb_{\star}, \tag{14}$ where $b_{\star} = \left[(G^\top \Sigma^{-1}G)C + I \right]^{-1}G^\top \Sigma^{-1}(y-Gm). \tag{15}$ It also holds that $C\left[(G^\top \Sigma^{-1}G)C + I \right]^{-1}G^\top \Sigma^{-1} = CG^\top \left(GCG^\top + \Sigma \right)^{-1}, \tag{16}$ which implies that (14) agrees with expression (6).

Basis Approach

We now consider an alternative approach that avoids the joint optimization over $(u,b)$. Similar exposition can be found in section 8.3 of (Evensen et al., 2022). Looking back to (9), instead of choosing to optimize over all $b$ satisfying (10), we will instead simply choose a particular $b$ satisfying this constraint. This only makes sense if the particular choice of $b$ does not matter; i.e., if the value of $\langle b, u^\prime \rangle$ is the same for any $b$ solving (10). The below propostion shows that this is indeed the case.

Proposition.
Let $u \in \mathbb{R}^d$ be a vector such that there is at least one solution $b \in \mathbb{R}^d$ to $Cb = u^\prime$. Then $\langle b, u^\prime \rangle$ and $J(u,b)$ are constant for any choice of $b$ solving this linear system.

Proof. Let $b, b^\prime \in \R^d$ satisfy $Cb = Cb^\prime = u^\prime$. It thus follows that $\langle b, u^\prime\rangle - \langle b^\prime, u^\prime\rangle = \langle b - b^\prime, u^\prime\rangle = \langle b - b^\prime, Cb\rangle = \langle C(b - b^\prime), b\rangle = \langle 0, b\rangle = 0,$ where we have used the linearity of the inner product and the fact that $C$ is symmetric. Since the inner product term in (11) is the only portion with dependence on $b$, it follows that $J(v,b) = J(v,b^\prime)$. $\qquad \blacksquare$

Thus, for each $u$ that yields a consistent system $Cb = u^\prime$, we can simply pick any solution $b$ to insert into $\langle b, u^\prime\rangle$. The objective will be well-defined since the above result verifies that the specific choice of $b$ is inconsequential. A natural choice is to select the $b$ of minimal norm; that is, $b^{\dagger} := \text{argmin}_{Cb=u^\prime} \lVert b \rVert. \tag{17}$ This unique minimal norm solution is guaranteed to exist and is conveniently given by the Moore-Penrose pseudoinverse $b^{\dagger} = C^{\dagger}u^\prime. \tag{18}$ Note that when $C$ is positive definite, $C^{\dagger} = C^{-1}$. We can now eliminate the requirement to optimize over $b$ by considering the following optimization problem.

Pseudoinverse Formulation.
\begin{align} u_{\star} &\in \text{argmin}_{u \in \mathcal{S}} J(u) \tag{19} \newline J(u) &:= \frac{1}{2} \lVert y - Gu\rVert^2_{\Sigma} + \frac{1}{2}\langle C^{\dagger}(u-m), u-m \rangle \newline \mathcal{S} &:= \left\{u \in \mathbb{R}^d: u-m \in \mathcal{R}(C) \right\} \end{align}

Note that we are now required to reintroduce a constraint set $\mathcal{S}$. This is due to the fact that the pseudoinverse produces a solution $b$ to $Cb = u^\prime$ when a solution exists, but we must still explicitly restrict the search space to $u$ that admit a consistent system $Cb = u^\prime$; i.e., $u^\prime \in \mathcal{R}(C)$. The previous approach dealt with this issue by introducing a Lagrange multiplier $\lambda$. We now show that (18) can be written as an unconstrained problem without extending the parameter space. This relies on the following important fact.

Subspace Property.
Let $A \in \mathbb{R}^{d \times r}$ be a matrix of rank $r \leq d$ satisfying $C = AA^\top$. Then a solution to the optimization problem (19) is constrained to the $r$-dimensional affine subspace $m + \mathcal{R}(C) = m + \mathcal{R}(A). \tag{20}$

This is simply a rewriting of the requirement $u^\prime = u-m \in \mathcal{R}(C)$. The solution is constrained by the rank of the subspace $\mathcal{R}(C)$. If $r = d$ (i.e., $C$ is full-rank) then we recover the typical unconstrained least squares problem. In general, $r < d$ so the columns of $C$ are linearly dependent. We therefore introduce the matrix $A$ whose columns $a_1, \dots, a_r$ form a basis for the range of $C$; i.e., the columns of $C$ and $A$ have the same span.

To implicitly encode the constraint $u^\prime \in \mathcal{R}(C)$, we can look for solutions in $m + \mathcal{R}(A)$. Any vector in this space can be written as $u = m + \sum_{j=1}^{r} w_j a_j = m + Aw, \tag{21}$ for some weight vector $w := (w_1, \dots, w_r)^\top \in \mathbb{R}^r$. We can now substitute $m + Aw$ for $u$ and $AA^\top$ for $C$ in (19) and thus formulate the optimization over the weights $w$: $J(w) := \frac{1}{2}\lVert y - G(m+Aw)\rVert^2_{\Sigma} + \frac{1}{2}\langle (AA^\top)^{\dagger}Aw, Aw \rangle. \tag{22}$ The optimization is now over $w \in \R^r$ so we have succeeded in writing (18) as an unconstrained problem. We can simplify this even further using the following properties of the pseudoinverse:

1. By definition, the pseudoinverse satisfies $A = A(A^{\dagger}A)$.
2. The matrix $A^{\dagger}A$ is a projection onto the rowspace of $A$; in particular it is idempotent (i.e., $(A^{\dagger}A)^2 = A^{\dagger}A$) and symmetric.
3. The following identity holds: $A^{\dagger} = A^\top (AA^\top)^{\dagger}$.

We can thus simplify the second term in (22) as \begin{align} \langle (AA^\top)^{\dagger}Aw, Aw\rangle = \langle A^\top(AA^\top)^{\dagger}Aw, w\rangle = \langle (A^{\dagger}A)w, w\rangle &= \langle (A^{\dagger}A)^2w, w\rangle \newline &= \langle (A^{\dagger}A)w, (A^{\dagger}A)w\rangle \newline &= \lVert (A^{\dagger}A)w \rVert^2, \tag{23} \end{align} where the second equality uses the third pseudinverse property above. The third and fourth equalities follow from the fact that $A^{\dagger}A$ is idempotent and symmetric. We can thus re-parameterize as $\tilde{w} := (A^{\dagger}A)w, \tag{24}$ and note that by the first pseudoinverse property we have $Aw = A(A^{\dagger}A)w = A\tilde{w}, \tag{25}$ which allows us to replace $Aw$ by $A\tilde{w}$ in the first term in (22). We obtain $J(\tilde{w}) := \frac{1}{2} \lVert y - G(m + A\tilde{w})\rVert_{\Sigma}^2 + \frac{1}{2} \lVert \tilde{w} \rVert^2 . \tag{26}$ We see that (26) is the sum of a model-data fit term plus a simple $L^2$ penalty on the weights. This final formulation is summarized below, where we have re-labelled $\tilde{w}$ as $w$ to lighten notation.

Unconstrained Pseudoinverse Formulation.
The optimization problem in (19) can equivalently be formulated as the unconstrained problem \begin{align} w_{\star} &\in \text{argmin}_{w \in \mathbb{R}^r} J(w) \tag{27} \newline J(w) &:= \frac{1}{2} \lVert y^\prime - GAw\rVert_{\Sigma}^2 + \frac{1}{2}\lVert w \rVert^2 \end{align}

One might initially take issue with the claim that (27) is actually unconstrained given that it relies on the re-parameterization (24), which constrains the weights to lie in the range of $A^{\dagger}A$. However, recalling that $A^{\dagger}A$ is an orthogonal projection, we have the following projection property: $\lVert \tilde{w} \rVert^2 = \lVert (A^{\dagger}A)w \rVert^2 \leq \lVert w \rVert^2. \tag{28}$ In other words, allowing the weights to be unconstrained can only increase the objective function (since switching $w$ and $\tilde{w}$ has no effect on the first term in (26)), so we are justified in considering the weights to be unconstrained in (27). In fact, we could have jumped right to this conclusion from (23) and avoided needing to define $\tilde{w}$ at all.

The following result provides the solution to (27), which immediately follows from the observation that (27) is a standard least squares problem.

Solution of Unconstrained Pseudoinverse formulation.
The minimizing weight of the optimization problem in (26) is given by $w_{\star} = A^\top G^\top \left(GCG^\top + \Sigma \right)^{-1}y^\prime, \tag{29}$ which implies the optimal parameter \begin{align} u_{\star} &= m + Aw_{\star} \newline &= m + CG^\top \left(GCG^\top + \Sigma \right)^{-1}(y - Gm) \end{align} agrees with the typical least squares solution in (6).

Proof. Observe that (27) is of the form (1) with $u := w$, $y^\prime := y$, $G := GA$, and $C := I$. Thus, we apply (6) to obtain \begin{align} w_{\star} &:= (GA)^\top \left[(GA)(GA)^\top + \Sigma \right]^{-1}y^\prime \newline &= A^\top G^\top \left[GCG^\top + \Sigma \right]^{-1}(y-Gm), \end{align} where the second equality uses $C = AA^\top$. We thus have \begin{align} u_{\star} &= m + Aw_{\star} \newline &= m + (AA^\top) G^\top \left[GCG^\top + \Sigma \right]^{-1}(y-Gm) \newline &= C G^\top \left[GCG^\top + \Sigma \right]^{-1}(y-Gm), \end{align} which exactly agrees with (6). $\qquad \blacksquare$

Concluding Notes

We have successfully shown that the optimum (6) continues to hold for the generalization (11) of the least squares formulation, even when $C$ is singular. We provided two different derivations of the solution, resulting in the expressions (14) and (29), both of which agree with $u_{\star}$ as given in (6). We conclude with a few brief follow-up comments.

Uniqueness

Our derivations show that the optimum $u_{\star}$ is unique, even when $C$ is singular. Note that this does not imply that a solution $(u_{\star}, b_{\star})$ to (11) is unique, since many choices of $b$ may lead to the same value of the objective. Indeed, $u_{\star} = m + Cb_{\star}$ but it may be that $u_{\star} = m + Cb^\prime$ for some other $b^\prime \in \mathcal{S}$. n particular, if $u_{\star}$ is optimal, then $(u_{\star}, b^\prime)$ minimizes $J(u,b)$ for any $b^\prime$ satisfying $Cb^\prime = u_{\star} - m$. This follows immediately from the result that the particular $b$ solving this linear system does not change the value of the objective. Given that the particular choice of $b$ is inconsequential, then a particular rule for choosing $b$ will not affect the optimal $u_{\star}$. The specific approach of choosing the minimal norm $b$ led to the least squares problem (27), which we know has a unique solution. Therefore, if (11) has a solution, then it is unique. We summarize this below.

Uniqueness.
If the solution set $\mathcal{S}$ is non-empty, then there is a unique optimal value $u_{\star}$ solving (11).

Representations of the solution

Continuing the above discussion, it should also be emphasized that, while the solution is unique, there may be many different ways to represent it. To see this, consider that (14) gives the solution in the form $u_{\star} = m + Cb_{\star}$, for a weight vector $b_{\star}$. But in general $C$ has linearly dependent columns, and thus there may be multiple sets of weights that give rise to the same vector. In (29) the solution is instead represented in the form $u_{\star} = m + Aw_{\star}$. If we assume that the columns of $A$ are linearly independent, thus providing a basis for $\mathcal{R}(C)$, then this provides the unique representation of the solution with respect to the particular basis $a_1, \dots, a_r$. Note that the two representations are connected by $u_{\star} = m + Cb_{\star} = m + AA^\top b_{\star} = m + A(A^\top b_{\star}) \tag{29}.$ If the columns of $A$ are independent, then (29) implies $w_{\star} = A^\top b_{\star}$.

Applications

Why should we even care about the case when $C$ is singular? The main application that motivated this post comes from the field of data assimilation. In this context, the parameter and data dimensions $d$ and $n$ can be massive, potentially precluding storing the matrices $\Sigma$ or $C$. A common solution to this problem is to replace the true $C$ with a low-rank approximation of the form $C = \frac{1}{r-1} \sum_{j=1}^{r} (u^{(j)} - m)(u^{(j)} - m)^\top, \tag{30}$ where $r < d$. This Monte Carlo approximation results in a sample covariance matrix $C$ that has rank at most $r$. Here, $\{u^{(j)}\}$ is a set of $r$ samples used to compute the sample covariance. If we define $A$ to be the $d \times r$ matrix with $j^{\text{th}}$ row equal to $\frac{1}{\sqrt{r-1}}\left(u^{(j)} - m\right)$, then we see that $C = AA^\top$. The matrix $A$ is full rank if the vectors $\{u^{(j)}\}$ are independent.

Appendix

Proof: Solution for Lagrange Multiplier Formulation

We derive expressions (14), (15) and (16), the latter showing that the solution agrees with (6).

Deriving (14) and (15)

We start by computing the gradients with respect to $u$, $b$, and $\lambda$: \begin{align} \nabla_{u}J &= -G^\top \Sigma^{-1}(y-Gu) + \frac{1}{2}b - \lambda \tag{31} \newline \nabla_{b}J &= \frac{1}{2}u^\prime + C\lambda \newline \nabla_{\lambda}J &= Cb - u^\prime. \end{align} where we recall the definition $u^\prime := u - m$. We now solve the system of equations $\nabla_{u}J = \nabla_{b}J = \nabla_{\lambda}J = 0$ for the three unknowns. Focusing on the last two equations, if we compute $2\nabla_{b}J + \nabla_{\lambda}J = 0$ we obtain the optimality criterion $C(b + 2\lambda) = 0, \tag{32}$ or equivalently, $b + 2\lambda \in \mathcal{N}(C), \tag{33}$ where $\mathcal{N}(C)$ denotes the null space of $C$. The null space of $C$ may be nontrivial, which aligns with the above discussion that many values of $b$ may yield the optimal solution. We consider taking the minimal norm solution $b_{\star} + 2\lambda_{\star} = 0, \tag{34}$ which implies $\lambda_{\star} = -\frac{1}{2}b_{\star}$. Note also that $y - Gu_{\star} = (y-Gm) - G(u_{\star}-m) =: y^\prime - Gu_{\star}^\prime = y^\prime - GCb_{\star}. \tag{35}$ We plug in the expression for $\lambda_{\star}$ to $\nabla_{u}J = 0$, which yields $-G^\top \Sigma^{-1}(y^\prime - GCb_{\star}) + b_{\star} = 0. \tag{36}$ Solvign for $b_{\star}$, we obtain $b_{\star} = \left[(G^\top \Sigma^{-1}G)C + I \right]^{-1} G^\top \Sigma^{-1}(y-Gm), \tag{35}$ which is equation (15). Recalling the definition $u^\prime_{\star} := u_{\star} - m$, we obtain equation (14): $u_{\star} = m + u^\prime_{\star} = m + Cb_{\star}. \tag{36}$

Deriving (16)

We now show that (36) agrees with the standard least squares solution (6). Differentiating the latter with a tilde, we want to show equality between the two expressions \begin{align} u_{\star} &= m + C\left[(G^\top \Sigma^{-1}G)C + I \right]^{-1} G^\top \Sigma^{-1}(y-Gm) \newline \tilde{u}_{\star} &= m + CG^\top \left(GCG^\top + \Sigma \right)^{-1}(y-Gm). \end{align} We therefore see that it suffices to show $\left[(G^\top \Sigma^{-1}G)C + I \right]^{-1} G^\top \Sigma^{-1} = G^\top \left[GCG^\top + \Sigma \right]^{-1}, \tag{37}$ which implies (16). To show (37), we multiply both sides by each of the matrices in brackets, which gives $G^\top \Sigma^{-1}\left[GCG^\top + \Sigma \right] = \left[G^\top \Sigma^{-1}GC + I \right] G^\top. \tag{38}$ The righthand side of (38) can be factored as $\left[G^\top \Sigma^{-1}GC + I \right] G^\top = G^\top \Sigma^{-1} \left[GCG^\top + \Sigma \right], \tag{39}$ completing the proof. $\qquad \blacksquare$