Generalized SVD for Multiple Penalties

Why Generalized Singular Coordinates Are Useful

For ordinary Tikhonov regularization, the penalty is proportional to the squared Euclidean norm of the parameter. When the penalty instead has the form \(\|Lz\|^2\), it can favor some directions in the parameter space and strongly penalize others. The singular vectors of the forward operator \(A\) alone may not reveal how those directions interact. The generalized singular value decomposition (GSVD) diagonalizes the data operator and penalty operator together, making the trade-off visible coordinate by coordinate.

We study the objective \(\|Az-y\|^2+\lambda\|Lz\|^2\), where \(y\) is the observed vector, \(z\) is the unknown parameter, \(A\) is the forward operator, \(L\) is the penalty operator, and \(\lambda>0\) controls the penalty strength. The essential coordinates are not generally the entries of \(z\) itself: they are the coordinates of \(z\) after a change of basis determined jointly by \(A\) and \(L\).

This normalization fixes a potential source of confusion. The generalized singular values \(\gamma_i\) are ratios, not substitutes for \(\alpha_i\) and \(\beta_i\) in every formula. In particular, the coordinatewise estimator multiplier must be derived using the normalized pair \((\alpha_i,\beta_i)\).

Existence and Coordinatewise Solution

The stated decomposition exists under the full-column-rank condition. One way to see this is to set \(B=A^\mathsf{T}A+L^\mathsf{T}L\), which is positive definite under that condition. Diagonalize the symmetric matrix \(B^{-1/2}A^\mathsf{T}AB^{-1/2}\) with an orthogonal matrix \(Q\), and set \(X=B^{-1/2}Q\). Then \(X^\mathsf{T}BX=I\), while \(X^\mathsf{T}A^\mathsf{T}AX\) is diagonal with entries \(\alpha_i^2\). Since \(X^\mathsf{T}L^\mathsf{T}LX=I-X^\mathsf{T}A^\mathsf{T}AX\), its diagonal entries are \(\beta_i^2=1-\alpha_i^2\). Orthonormal columns \(U,V\) can be chosen to represent \(AX\) and \(LX\), including by completing the columns when a diagonal entry is zero.

Proof. Write \(w=X^{-1}z\). Then \(z=Xw\), and the two operator representations give \(Az=UCw\) and \(Lz=VSw\). Because the columns of \(U\) are orthonormal, \(\|UCw-y\|^2=\|Cw-U^\mathsf{T}y\|^2+\|y-UU^\mathsf{T}y\|^2\). The second term does not depend on \(w\). Orthonormality of \(V\) also gives \(\|VSw\|^2=\|Sw\|^2\). Therefore, apart from the constant residual term, the objective is \(\sum_{i=1}^n[(\alpha_iw_i-u_i^\mathsf{T}y)^2+\lambda\beta_i^2w_i^2]\). For each \(i\), differentiating this scalar quadratic with respect to \(w_i\) gives \(2\alpha_i(\alpha_iw_i-u_i^\mathsf{T}y)+2\lambda\beta_i^2w_i=0\). Thus \((\alpha_i^2+\lambda\beta_i^2)w_i=\alpha_i u_i^\mathsf{T}y\). The denominator is positive because \(\alpha_i^2+\beta_i^2=1\) and \(\lambda>0\). Each scalar quadratic is strictly convex, so these coordinates give the unique minimizer. Finally, when \(\beta_i>0\), normalization implies \(\alpha_i=\gamma_i/\sqrt{1+\gamma_i^2}\) and \(\beta_i=1/\sqrt{1+\gamma_i^2}\). Substitution yields \(\dfrac{\alpha_i}{\alpha_i^2+\lambda\beta_i^2}=\dfrac{\gamma_i\sqrt{1+\gamma_i^2}}{\gamma_i^2+\lambda}\). This proves the result. \(\square\)

The fitted values have an especially direct form. Multiplying the estimated coordinates by \(A\) gives \(A\hat z=U\operatorname{diag}\!\left(\dfrac{\alpha_i^2}{\alpha_i^2+\lambda\beta_i^2}\right)U^\mathsf{T}y\). Thus, in the data-space directions \(u_i\), the fitted-value multipliers are \(\alpha_i^2/(\alpha_i^2+\lambda\beta_i^2)\). This filter is different from the multiplier for estimating \(w_i\): the latter includes one factor \(\alpha_i\), rather than \(\alpha_i^2\), in its numerator.

What the Generalized Singular Values Say

When \(\beta_i>0\), the coordinate multiplier can be expressed using \(\gamma_i\), but the normalization factors matter: \(\alpha_i/(\alpha_i^2+\lambda\beta_i^2)=\gamma_i\sqrt{1+\gamma_i^2}/(\gamma_i^2+\lambda)\). The fitted-value multiplier has the simpler ratio form \(\alpha_i^2/(\alpha_i^2+\lambda\beta_i^2)=\gamma_i^2/(\gamma_i^2+\lambda)\). A large \(\gamma_i\) means the data operator is strong relative to the penalty in that generalized direction, so its fitted-value multiplier is closer to one. A small \(\gamma_i\) means the penalty has greater relative influence.

There are also endpoint cases. If \(\beta_i=0\), then \(\alpha_i=1\): that direction is unpenalized, and its transformed estimate is \(u_i^\mathsf{T}y\). The ratio \(\gamma_i\) is not finite in this case, so use the \((\alpha_i,\beta_i)\) formula. If \(\alpha_i=0\), then \(\beta_i=1\): the data do not identify that direction, and the minimizing transformed coordinate is zero. The full-column-rank condition on the stacked operator rules out \(\alpha_i=\beta_i=0\), which would leave that parameter direction unidentified by both operators.

Several Penalties and the Combined Operator

A problem with several quadratic penalties can be put into the same form. Suppose \(L_1,\ldots,L_k\) are penalty operators and \(\lambda_j>0\). Define the vertically stacked, weighted operator \(P=\begin{bmatrix}\sqrt{\lambda_1}L_1\\ \vdots\\ \sqrt{\lambda_k}L_k\end{bmatrix}\). Then \(\|Pz\|^2=\sum_{j=1}^k\lambda_j\|L_jz\|^2\), so minimizing \(\|Az-y\|^2+\sum_j\lambda_j\|L_jz\|^2\) is exactly a two-operator Tikhonov problem with penalty operator \(P\) and penalty coefficient one. If the stacked operator \(\begin{bmatrix}A\\P\end{bmatrix}\) has full column rank and the dimensions permit the stated GSVD, the theorem applies to \(A\) and \(P\).

Interpretation and a Common Mistake

The GSVD is useful because it separates the interaction between fit and penalty into scalar trade-offs. It is especially important when the penalty is not proportional to the identity: a standard SVD of \(A\) alone generally does not diagonalize \(\|Lz\|^2\) in the same coordinates. The GSVD instead supplies a single transformed coordinate system in which both quadratic terms are sums of independent coordinate contributions.

Common mistake: Starting from \(\alpha_i^2+\beta_i^2=1\) and \(\gamma_i=\alpha_i/\beta_i\), it is tempting to replace \(\alpha_i/(\alpha_i^2+\lambda\beta_i^2)\) by \(\gamma_i/(\gamma_i^2+\lambda)\). But normalization gives \(\alpha_i=\gamma_i/\sqrt{1+\gamma_i^2}\) and \(\beta_i=1/\sqrt{1+\gamma_i^2}\). Substitution therefore yields \(\gamma_i\sqrt{1+\gamma_i^2}/(\gamma_i^2+\lambda)\), not \(\gamma_i/(\gamma_i^2+\lambda)\). The latter expression omits a factor and fails, for instance, at \(\gamma_i=\lambda=1\). A reliable check is to derive the multiplier first in \(\alpha_i,\beta_i\), then substitute the ratio and verify a simple numerical case.

A second distinction worth keeping explicit is that \(\hat w_i\) is a coordinate in the transformed parameterization \(z=Xw\), not usually the \(i\)th entry of the original parameter vector. The original estimate is \(X\hat w\). Confusing these vectors can yield an apparently plausible but incorrectly scaled answer.