How to compute the "hat-matrix" of constrained least squares

Question

I'm attempting to calculate the studentized residuals on a (equality) constrained least-squares regression for outlier detection. However, i'm a little uncertain on how to calculate the leverages, $h_i$, which is the diagonal of the "hat-matrix", $H$.

In the ordinary least squares case (without linear constraints), this matrix is computed as

$$ H = X\left(X^T X\right)^{-1}X^T \quad \quad \quad (1) $$

Where $X$ is the design matrix and $y$ is the variable to be explained. The constrained least-squares regression i'm attempting to do is the following:

$$ \left[ \begin{array}{c} \hat{\beta} \\ \hat{\lambda} \end{array} \right] = \left[ \begin{array}{cc} 2 X^T X & C^T \\ C & 0 \end{array} \right]^{-1} \left[ \begin{array}{c} 2 X^T y \\ d \end{array} \right] $$

where $X$ is again the design matrix, $y$ the variable to be explained and $C$ is the restriction matrix such that

$$ C \beta = d. $$

$\lambda$ is in this case the lagrange multipliers. My Question is then, what is the "hat"-equivalent matrix for this type of regression? Will the formulation in $(1)$ hold? My quess is no, since you are adding additional information to the regression. For instance, if you added 0-restrictions to some of the columns, you might as well have excluded them from the design matrix, in which case the leverages would change.

Answer 1

It's been quite a while since this question was asked, but searching the internet still gives only few information.

For unconstrained weighted linear least squares regression (WLS) with the design matrix $X$ ($n$ x $m$), the response vector $y$ ($n$ x $1$) and the weight-matrix $W = diag\left(w_i\right)$ with the weight-vector $w$ ($n$ x $1$) assigning one weight to each sample $\left(x_i^T\ |\ y_i\right)$, the "hat" or projection matrix $H$ ($n$ x $n$) is defined as

$$ H = X\left(X^{T}WX\right)^{-1}X^{T}W \quad\quad\quad\quad\quad\left(1\right) $$

and it is called projection matrix because it projects the response vector $y$ to obtain its respective prediction vector $\hat{y}$ by

$$ \hat{y} = Hy \quad\quad\quad\quad\quad\left(2\right) $$

which could also be obtained by

$$ \hat{y} = Hy = X\left(X^{T}WX\right)^{-1}X^{T}Wy = X\hat{\beta} \quad\quad\quad\quad\quad\left(3\right) $$

using the coefficient vector $\hat{\beta}$ ($m$ x $1$)

$$ \hat{\beta} = \left(X^{T}WX\right)^{-1}X^{T}Wy \quad\quad\quad\quad\quad\left(4\right) $$

(Note that in the original question $y$ is included at the end of $H$ whereas $X$ should be included at the beginning). For the unconstrained ordinary least squares (OLS) problem You have stated, the weight-matrix $W$ is simply the identity matrix $I$.

The definition $\left(1\right)$ is no longer valid as soon as the constraints are introduced, but $\left(2\right)$ is still applicable as it only states the general projection operation, which one can examine in more detail. According to $\left(2\right)$ $\hat{y}$ is a linear combination of the responses in $y$ with each row in $H$ holding the respective weights for these linear combinations

$$ \hat{y}_{i} = \sum_{j=1}^{n}{\left(h_{ij}\ \cdot\ y_{j}\right)} = h_{i}^{T} y \quad\quad\quad\quad\quad\left(5\right) $$

Therefore, one could also formulate that the entry $h_{ij}$ of the hat-matrix (which is $= h_{ji}$ only for OLS) is the slope that determines how much $\hat{y}_{i}$ changes as $y_{j}$ is altered. In other words $h_{ij}$ can be interpreted as a derivative

$$ h_{ij} = \frac{\partial \hat{y}_{i}}{\partial y_{j}} \quad\quad\quad\quad\quad\left(6\right) $$

and consequently

$$ H = \frac{\partial \hat{y}}{\partial y} \quad\quad\quad\quad\quad\left(7\right) $$

Combining $\left(7\right)$ and $\left(3\right)$, this can also be re-expressed as

$$ \frac{\partial \hat{y}}{\partial y} = H = \frac{\partial X\hat{\beta}}{\partial y} = \frac{\partial X}{\partial y}\hat{\beta} + X\frac{\partial \hat{\beta}}{\partial y} = 0\hat{\beta} + X\frac{\partial \hat{\beta}}{\partial y} = X\frac{\partial \hat{\beta}}{\partial y} \quad\quad\quad\quad\quad\left(8\right) $$

(see the Matrix cookbook chapter 2 (37)).
Again, this only relies on the general projection operation and is not a special case derived for WLS or OLS. Now, as You have written, the constraint $C\hat{\beta} = d$ ($C$ has shape $o$ x $m$) can be included, the Lagrangian gradient can be set to 0 and the resulting system of equations is solved by

$$ \left[ \begin{array}{c} \hat{\beta} \\ \hat{\lambda} \end{array} \right] = \left[ \begin{array}{cc} X^T W X & C^T \\ C & 0 \end{array} \right]^{-1} \left[ \begin{array}{c} X^T W y \\ d \end{array} \right] = Q^{-1} \left[ \begin{array}{c} X^T W y \\ d \end{array} \right] \quad\quad\quad\quad\quad\left(9\right) $$

(I've multiplied the least squares objective with 0.5 to get rid off the factor 2, which does not affect the result)
If one denotes

$$ \left[ \begin{array}{c} \hat{\beta} \\ \hat{\lambda} \end{array} \right] $$

by $\beta^{\ast}$ and introduces the augmented matrix ($n$ x ($m$ + $o$))

$$ X^{\ast} = \left[ \begin{array}{cc} X & 0 \end{array} \right] $$

where the last $o$ columns are $0$-vectors for the Lagrangian multipliers in $\hat{\lambda}$, predictions can be made by

$$ \hat{y} = X^{\ast}\beta^{\ast} \quad\quad\quad\quad\quad\left(10\right) $$

Applying $\left(8\right)$ results in

$$ H = \frac{\partial \hat{y}}{\partial y} = X^{\ast}\frac{\partial \beta^{\ast}}{\partial y} \quad\quad\quad\quad\quad\left(11\right) $$

So, $\left(9\right)$ needs to be derived with respect to $y$ which gives

$$ \frac{\partial \beta^{\ast}}{\partial y} = Q^{-1} \left[ \begin{array}{c} X^T W \\ 0 \end{array} \right] \quad\quad\quad\quad\quad\left(12\right) $$

(assuming $d \neq f\left(y\right)$; see the Matrix cookbook chapter 2 (37)).
This can be continued to give the final solution

$$ H = X^{\ast} Q^{-1} \left[ \begin{array}{c} X^T W \\ 0 \end{array} \right] \quad\quad\quad\quad\quad\left(13\right) $$

For computation, it might be more convenient to form $Q^{-1\ast}$ which is simply $Q^{-1}$ with the last $o$ rows and columns skipped ($Q^{-1}$ has shape (($m$ + $o$) x ($m$ + $o$))). This does not imply that the corresponding entries need not to be calculated. Since their precence affects $Q^{-1}$ overall, they may only be skipped to ommit the matrix multiplication of $X^{\ast}$ with the Lagrangian multipliers $\hat{\lambda}$ in $\beta^{\ast}$.
With this, the simplified formulation is

$$ H = X Q^{-1\ast} X^{T}W \quad\quad\quad\quad\quad\left(14\right) $$

Comparing $\left(14\right)$ and $\left(1\right)$, the only difference that can be observed is the change of the inverse from $\left(X^{T}WX\right)^{-1}$ to $Q^{-1\ast}$, so $\left(14\right)$ seems to be a generalization of $\left(1\right)$ and indeed if the equality constraint $C\hat{\beta} = d$ is removed, $Q^{-1\ast}$ collapses to $\left(X^{T}WX\right)^{-1}$ as expected for the unconstrained case.

I have tested this for an equality constrained ridge regression for which the Lagrangian system is solved by

$$ \left[ \begin{array}{c} \hat{\beta} \\ \hat{\lambda} \end{array} \right] = \left[ \begin{array}{cc} X^T W X + \alpha I& C^T \\ C & 0 \end{array} \right]^{-1} \left[ \begin{array}{c} X^T W y \\ d \end{array} \right] = Q^{-1} \left[ \begin{array}{c} X^T W y \\ d \end{array} \right] $$

with the regularization parameter $\alpha$ that allows to set a certain sum of squared residuals for the fit (smoothing; basically just another generalization). For this purpose, I have compared if

$$ Hy = X\hat{\beta} $$

holds (both give $\hat{y}$) and they have coincided well for the tests I conducted. Besides, I had $h_{ii} \leq 1$.

For WLS, $H$ is only symmetric if all weights are equal (= OLS). Otherwise, this is not the case because a response $y_i$ with a small weight will have little effect on the response $y_j$ with a high weight ($h_{ji}$ is small) whereas $y_j$ will exhibit a large influence on $y_i$ ($h_{ij}$ is large and thus $\neq h_{ji}$). Also, the idempotence

$$ HH = H $$

will no longer be given for the constrained case in $\left(14\right)$. Without the constraints, it would be

$$ HH = X\left(X^{T}WX\right)^{-1}X^{T}W X\left(X^{T}WX\right)^{-1}X^{T}W = X\left(X^{T}WX\right)^{-1}X^{T}W = H $$

but as $Q^{-1\ast}$ is not necessarily the inverse of $X^{T}WX$ the constrained least squares fit will act as a smoother with subsequent applications of $H$ as $HHH\cdots Hy$ always smoothing the prediction $\hat{y}$ a little further because the constrained least squares solution is not the unconstrained one.

Regarding the validity range of $H$ in $\left(14\right)$:
I have no method at hand to check the validity range of $H$, but I took a look at

Riazoshams, et al. "On the outlier Detection in Nonlinear Regression", International Journal of Mathematical and Computational Sciences Vol:3, No:12, 2009

where the hat-matrix is constructed using the slopes of the tangent hyperplane to replace $X$ in $\left(1\right)$ and depending of the curvature of the nonlinear regression problem, this can have a very limited validity range, too. So, for a basic outlier treatment, $\left(14\right)$ should be sufficient. I would recommend to work with Cook's Distance and/or Peña's Sensitivity from here.

Unfortunately, I'm not used to SAS. Thus I cannot check if it uses a similar definition or not.

Edit
Since the partial derivative of $\left(14\right)$ with respect to $\hat{\beta}$ is zero (since $\hat{\beta}$ is not part of the equation), this formulation of the hat matrix has no limit in validity. In other words, it is not a local function dependig on $\hat{\beta}$ as encountered in nonlinear regression, but applies globally as in Ordinary Linear Least Squares.