Tactics for Improving Least Squares Estimation

Here are the notational conventions used throughout this article. All vectors and matrices appear in boldface. An entry of a vector $\bm{x}$ or matrix $\bm{A}$ is denoted by the corresponding subscripted lower-case letter $x_{i}$ or $a_{ij}$. All entries of the vector $\bm{0}$ equal $0$; $\bm{I}$ indicates an identity matrix. The ^⊤ superscript indicates a vector transpose. The symbols $\bm{x}_{m}$ and $\bm{A}_{m}$ denote a sequence of vectors and a sequence of matrices with entries $x_{mi}$ and $a_{mij}$. The Euclidean norm and $\ell_{1}$ norm of a vector $\bm{x}$ are denoted by $\|\bm{x}\|$ and $\|\bm{x}\|_{1}$. The spectral, Frobenius, and nuclear norms of a matrix $\bm{A}$ are denoted by $\|\bm{A}\|$, $\|\bm{A}\|_{F}$, and $\|\bm{A}\|_{*}$. For a smooth real-valued function $f(\bm{\beta})$, we write its first differential (row vector of partial derivatives) as $df(\bm{\beta})$. The gradient $\nabla f(\bm{\beta})$ is the transpose of $df(\bm{\beta})$. The directional derivative of $f(\bm{\beta})$ in the direction $\bm{v}$ is written $d_{\bm{v}}f(\bm{\beta})$. When the function $f(\bm{\beta})$ is differentiable, $d_{\bm{v}}f(\bm{\beta})=df(\bm{\beta})\bm{v}$. The second differential (Hessian matrix) of $f(\bm{\beta})$ is $d^{2}f(\bm{\beta})$.

In M estimation (de Menezes et al., 2021; Huber and Dutter, 1974; Huber and Ronchetti, 2011) based on residuals, one minimizes a function

of the residuals $r_{i}=y_{i}-\bm{x}_{i}^{\top}\bm{\beta}$ with $\rho(r)$ bounded below. We have already dealt with the choice $\rho(r)=|r|$ in LAD regression. In general, one can replace $\rho(r)$ by its Moreau envelope, leverage the Moreau majorization, and immediately invoke ordinary least squares. Many authors advocate estimating $\bm{\beta}$ by minimizing a function

of the squared residuals. Conveniently, this functional form is consistent with loglikelihoods under an elliptically symmetric distribution (Lange and Sinsheimer, 1993). If $\rho(s)$ is increasing, differentiable, and concave, then

for all $s$. This plays out as the majorization

with weights $w_{mi}=\rho^{\prime}[(y_{i}-\bm{x}_{i}^{\top}\bm{\beta}_{m})^{2}]$ and an irrelevant constant $c_{m}$ that depends on $\bm{\beta}_{m}$ but not on $\bm{\beta}$. This surrogate can be minimized by deweighting and invoking ordinary least squares. The Geman-McClure choice $\rho(s)=\log(1+s)$ and the Cauchy choice $\rho(s)=\frac{s}{1+s}$ for $s\geq 0$ show that the intersection of the imposed conditions is non-empty.

A prime example (Chi and Chi, 2022; Heng et al., 2023) of robust M estimation revolves around the $\text{L}_{2}\text{E}$ objective (Scott, 2001; Liu et al., 2023)

This loss is particularly attractive because it incorporates a precision parameter $\tau$ as well as the regression coefficient vector $\bm{\beta}$. Because $-e^{-s}$ is concave, one can construct the MM surrogate

at iteration $m$ with $\tau$ fixed. The resulting weighted least squares surrogate can be further deweighted as explained in Section 2.5. Under deweighting it morphs into the surrogate

where $w_{mi}=e^{-\tau^{2}r_{mi}^{2}/2}$ and $c_{m}$ is an irrelevant constant. This ordinary least squares surrogate allows one to recycle a single Cholesky decomposition across all iterations. For $\bm{\beta}$ fixed, one can update the precision parameter $\tau$ by gradient descent (Heng et al., 2023) or an approximate Newton’s method (Liu et al., 2023). Both methods rely on backtracking.

Quantile regression (Koenker and Bassett, 1978) is a special case of M estimation with the objective $f(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}\rho_{q}(y_{i}-\bm{x}_{i}^{\top}\bm{\beta})$, where

is called the check function. The argument $r$ here suggests a residual. When $q=\frac{1}{2}$, quantile regression reduces to LAD regression. The check function is not differentiable at $r=0$. To smooth the objective (9), one can take either its Moreau envelope or the Moreau envelope of just the nonsmooth part $\frac{1}{2}|r|$. In the former case, elementary calculus shows that the minimum is attained at

It follows that

A benefit of replacing $\rho_{q}(r)$ by $M_{\mu\rho_{q}}(r)$ is that, by definition, $M_{\mu\rho_{q}}(r)$ satisfies the Moreau majorization (4) at the anchor point $r_{m}$. This majorization can be improved, but it is convenient because it induces the uniform majorization

which is an unweighted sum of squares in the residuals $r_{i}=y_{i}-\bm{x}_{i}^{\top}\bm{\beta}$.

Convolution approximation offers yet another avenue for majorization (Fernandes et al., 2021; He et al., 2023; Kaplan and Sun, 2017; Tan et al., 2022; Whang, 2006). If we select a well-behaved kernel $k(x)\geq 0$ with bounded support and total mass $\int k(x)dx=1$, then the convolution

approximates $g(y)$ as the bandwidth $\mu\downarrow 0$. Given $g(x)=|x|$ and the symmetric uniform kernel $k(x)=\frac{1}{2}1_{\{|x|\leq 1\}}$, the necessary integral

is straightforward to calculate. This is just the Moreau envelope $M_{\mu|\cdot|}(y)$ shifted upward by $\frac{\mu}{2}$. This approximate loss is now optimized via the Moreau majorization

where $c_{m}$ is an irrelevant constant and $z_{m}=\mathop{\rm prox}\nolimits_{\mu|\cdot|}(r_{m})$. Although our numerical experiments in quantile estimation feature the convolution-smoothed quantile loss (11), the choice (10) works equally well. What is more important in practice is the transformation of quantile regression into a sequence of ordinary least squares problems with shifted responses.

The $\ell_{0}$-norm, the rank function, and their corresponding constraint sets are helpful in inducing sparsity and low matrix rank. Directly imposing these nonconvex penalties or constraint sets can lead to intractable optimization problems. A common middle ground is to use convex surrogates, such as the $\ell_{1}$-norm and the nuclear norm, but these induce shrinkage bias and a surplus of false positives in support recovery. The Moreau envelopes of these nonconvex penalties offer an attractive alternative that leverages differentiability and majorization.

The $\ell_{0}$ norm has Moreau envelope (Beck, 2017)

The corresponding proximal map sends $\beta_{j}$ to itself when $\frac{1}{2}\beta_{j}^{2}\geq\mu$ and to $0$ when $\frac{1}{2}\beta_{j}^{2}<\mu$. The Moreau envelope of the $0/\infty$ indicator of the sparsity set $S_{k}=\{\bm{\beta}:\|\bm{\beta}\|_{0}\leq k\}$ is determined by projection onto $S_{k}$. Elementary arguments show that

where $\beta_{(1)},\ldots,\beta_{(p-k)}$ denote the $p-k$ smallest $\beta_{j}$ in magnitude (Beck, 2017). The Moreau majorization defined by equation (4) is available in both examples. In general, the Moreau envelope of the $0/\infty$ indicator of a set $S$ is $\frac{1}{2\mu}\mathop{\rm dist}\nolimits(\bm{y},S)^{2}$, where $\mathop{\rm dist}\nolimits(\bm{y},S)$ is the Euclidean distance from $\bm{y}$ to $S$. Other nonconvex sparsity inducing penalties $p(\beta)$, such as the SCAD penalty (Fan and Li, 2001) and the MCP penalty (Zhang, 2010), also have straightforward Moreau envelopes (Polson et al., 2015).

For matrices, low-rank is often preferred to sparsity in estimation and inference. Let $R_{k}$ denote the set $p\times q$ matrices of rank $k$ or less. The Moreau envelope of the $0/\infty$ indicator of $R_{k}$ turns out to be $\frac{1}{2\mu}\mathop{\rm dist}\nolimits(\bm{B},R_{k})^{2}=\frac{1}{2\mu}\sum_{j>k}\sigma_{j}^{2}$, where the $\sigma_{j}$ are the singular values of $\bm{B}$ in descending order. If $\bm{B}$ has singular value decomposition $\bm{U}\bm{\Sigma}\bm{V}^{\top}$, then the corresponding proximal map takes $\bm{B}$ to $\bm{U}\bm{\Sigma}_{k}\bm{V}^{\top}$, where $\bm{\Sigma}_{k}$ sets the singular values of $\bm{\Sigma}$ smaller than $\sigma_{k}$ to $0$. This result is just the content of the Eckart-Young theorem. The function $\mathrm{rank}(\bm{B})$ has Moreau envelope (Hiriart-Urruty and Le, 2013)

The map $\mathop{\rm prox}\nolimits_{\mu\,\mathrm{rank}}(\bm{B})=\bm{U}\bm{\Sigma}_{k}\bm{V}^{\top}$ retains $\sigma_{j}$ when $\frac{1}{2}\sigma_{j}^{2}\geq\mu$ and sets it $0$ otherwise. In other words, the proximal operator hard thresholds the singular values. In both cases the Moreau majorization (4) is available.

The nuclear norm $\|\bm{B}\|_{*}=\sum_{i}\sigma_{i}$ (Recht et al., 2010) is a widely adopted convex penalty in low-rank estimation. Its proximal map soft thresholds all $\sigma_{i}$ by subtracting $\mu$ from each, then replacing negative values by $0$, and finally setting $\mathop{\rm prox}\nolimits_{\mu\|\cdot\|_{*}}(\bm{B})=\bm{U}\tilde{\bm{\Sigma}}\bm{V}^{\top}$, where $\tilde{\bm{\Sigma}}$ contains the thresholded values. The Moreau envelope $M_{\mu\|\cdot\|_{*}}(\bm{B})$ serves as a smooth-approximation to $\|\bm{B}\|_{*}$ and benefits from the Moreau majorization (4).

To remove the weights and shift the responses in least squares, Heiser (1987) and Kiers (1997) suggest a crucial deweighting majorization. At iteration $m$, assume that the weights satisfy $0\leq w_{mi}\leq 1$ and that $\mu_{i}$ abbreviates the mean $\mu_{i}(\bm{\beta})$ of case $i$. Then we have

where $d_{m}=\sum_{i=1}^{n}w_{mi}(1-w_{mi})(y_{i}-\mu_{mi})^{2}$ does not depend on the current $\mu_{i}$. Equality holds, as it should, when $\mu_{i}=\mu_{mi}$.

The quadratic upper bound principle applies to functions $f(\bm{\beta})$ with bounded curvature (Böhning and Lindsay, 1988). Given that $f(\bm{\beta})$ is twice differentiable, we seek a matrix $\bm{H}$ satisfying $\bm{H}\succeq d^{2}f(\bm{\beta})$ and $\bm{H}\succ{\bf 0}$ in the sense that $\bm{H}-d^{2}f(\bm{\beta})$ is positive semidefinite for all $\bm{\beta}$ and $\bm{H}$ is positive definite. The quadratic bound principle then amounts to the majorization

The quadratic surrogate is exactly minimized by the update (2) with $g(\bm{\beta}_{m}\mid\bm{\beta}_{m})$ replaced by $\bm{H}$. The quadratic lower bound principle involves minorization and subsequent maximization.

The quadratic upper bound principle applies to logistic regression and its generalization to multinomial regression (Böhning, 1992; Krishnapuram et al., 2005). In multinomial regression, items are draw from $c$ categories numbered $1$ through $c$. Logistic regression is the special case $c=2$. If $y_{i}$ denotes the response of sample $i$, and $\bm{x}_{i}$ denotes a corresponding predictor vector, then the probability $w_{ij}$ that $y_{i}=j$ is

where $r_{ij}=\bm{x}_{i}^{\top}\bm{\beta}_{j}$. To estimate the regression coefficient vectors $\bm{\beta}_{j}$ over $n$ independent trials we must find a quadratic lower bound for the loglikelihood

Here $\bm{B}$ is the matrix with $j$th column $\bm{\beta}_{j}$. It suffices to find a quadratic lower bound for each sample, so we temporarily drop the subscript $i$ to simplify notation.

The loglikelihood has directional derivative

for the perturbation $\bm{V}$ of $\bm{B}$ with $j$th column $\bm{v}_{j}$. In the Kronecker product representation of $d_{\bm{V}}\mathcal{L}(\bm{B})$, $\bm{y}$ is now an indicator vector of length $c-1$ with a $1$ in entry $j$ if the sample belongs to category $j$, and $\bm{w}$ is the weight vector with value $w_{j}=\frac{e^{r_{j}}}{1+\sum_{k=1}^{c-1}e^{r_{k}}}$ in entry $j<c$. The quadratic form associated with the second differential $d^{2}\mathcal{L}(\bm{B})$ is

Across all cases Böhning (1992) derives the lower bound

where the $(c-1)\times(c-1)$ matrix $\bm{E}=\frac{1}{2}(\bm{I}-\frac{1}{c}\bm{1}\bm{1}^{\top})$ has the explicit inverse $2(\bm{I}+\bm{1}\bm{1}^{\top})$. The full-sample multinomial loglikelihood is then minorized by the quadratic

Maximizing the surrogate yields the update

Reverting to matrices, this update can be restated as

where $\bm{W}_{m}\in\mathbb{R}^{n\times(c-1)}$ conveys the weights for the first $c-1$ categories and the rows of $\bm{Y}\in\mathbb{R}^{n\times(c-1)}$ are indicator vectors conveying category membership. In the special case of logistic regression, the solution reduces to

Later, when we discuss penalized multinomial regression, we will consider the penalized loglikelihood $\mathcal{L}(\bm{B})-\frac{\lambda}{2}\|\bm{B}-\bm{P}_{m}\|_{F}^{2}$, where the projection $\bm{P}_{m}$ depends on the current iterate $\bm{B}_{m}$ and $\lambda>0$. Redoing the minorization (13) with this extra term leads to the stationary condition

with $\bm{C}=\bm{X}^{\top}(\bm{Y}-\bm{W}_{m})+\lambda\bm{P}_{m}$ and $\bm{\Delta}=\bm{B}-\bm{B}_{m}$. This condition is the same as the two equivalent matrix equations