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

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.

• The diagonal of $\hat H$ provides leverages only when constraints do not hold. In other cases you need to begin with an adequate definition of "leverage." Although many articles confuse the two (such as the Wikipedia article on leverage!), they are not the same. See en.wikipedia.org/wiki/… for a better definition based on influence: "Leverage is a measure of how much the estimated value of the dependent variable changes when the point is removed."
– whuber
Nov 14, 2018 at 15:36
• Right. That makes sense. I'm benchmarking my calculation with one made in SAS, and i'm quite close in getting the same "leverages" as it produces, but not quite there. I can't find any information on how SAS calculates these under constraints, though.
– amri
Nov 14, 2018 at 15:54
• This answer (see paragraphs around Eq. (4)) should be helpful. Apr 30 at 1:31

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.