Estimating SAVE With Conditional Covariance Matrices

From Conditional Covariances to a Sample Operator

SAVE uses how the predictor covariance changes across response values. Its sample version therefore needs more than slice means: within each slice, it estimates a covariance matrix, compares that matrix with the standardized overall covariance, and aggregates the squared differences. The resulting matrix is symmetric, so its eigenvectors provide the estimated directions.

Let $x_i\in\mathbb{R}^p$ be the predictor vector and $y_i$ the response for observation $i$, for $i=1,\ldots,n$. Partition the response values into $H$ nonempty slices. Write $\mathcal{I}_h$ for the indices in slice $h$, and $n_h$ for its count. This tutorial uses the empirical covariance within a slice with divisor $n_h$, and requires $n_h\geq 2$ for every slice.

Definition: The sample SAVE operator is formed after standardizing the predictors to have overall sample covariance equal to the identity matrix. If $C_h$ is the empirical covariance matrix of the standardized predictors in slice $h$, and $p_h=n_h/n$, then $$ \widehat{\mathrm{SAVE}}=\sum_{h=1}^{H}p_h(I_p-C_h)^2. $$ Here $I_p$ is the $p\times p$ identity matrix, and the square is matrix multiplication: $(I_p-C_h)^2=(I_p-C_h)(I_p-C_h)$.

Standardize Before Estimating Slice Covariances

Compute the sample mean and the usual sample covariance matrix of the predictors:

\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i, \qquad S=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^T.

Assume $S$ is positive definite. Its symmetric inverse square root $S^{-1/2}$ is defined by taking the inverse square roots of its positive eigenvalues. Standardize each observation by

z_i=S^{-1/2}(x_i-\bar{x}).

The standardized sample has mean zero and sample covariance $I_p$. Indeed, its sample mean is $S^{-1/2}(\bar{x}-\bar{x})=0$. Its sample covariance is

\begin{aligned} \frac{1}{n-1}\sum_{i=1}^n z_i z_i^T &=S^{-1/2} \left\{\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^T\right\} S^{-1/2}\\ &=S^{-1/2}SS^{-1/2}=I_p. \end{aligned}

This identity explains why SAVE compares each slice covariance with $I_p$. If $S$ is singular, this standardization is not defined in all $p$ directions; a lower-dimensional representation or a separately specified regularization is then needed.

For each slice, calculate its standardized mean and covariance:

\bar{z}_h=\frac{1}{n_h}\sum_{i\in\mathcal{I}_h}z_i, \qquad C_h=\frac{1}{n_h}\sum_{i\in\mathcal{I}_h} (z_i-\bar{z}_h)(z_i-\bar{z}_h)^T.

Each $C_h$ is symmetric and positive semidefinite. The weight $p_h=n_h/n$ is the observed fraction of the sample in slice $h$; the weights are nonnegative and sum to one. The divisor $n_h$ treats a slice as its empirical conditional distribution. Replacing it with $n_h-1$ gives a different finite-sample convention, so a calculation must state which convention it uses.

Worked Example: Standardizing and Computing Two Slice Variances

Suppose the one-dimensional predictor observations are $0,2,4,6$, and the first and last observations form one response slice while the middle two form the other. Standardize using the sample covariance divisor $n-1$, then calculate the sample SAVE estimate using divisor $n_h$ within each slice.

The sample mean is $\bar{x}=3$, and $S=\{(-3)^2+(-1)^2+1^2+3^2\}/3=20/3$. Let $a=\sqrt{3/20}$. The standardized observations, in order, are $-3a,-a,a,3a$.

Each slice has weight $1/2$. The first slice has standardized values $-3a,3a$, mean zero, and covariance $C_1=\{(3a)^2+(3a)^2\}/2=9a^2=27/20$. The second slice has values $-a,a$, mean zero, and covariance $C_2=(a^2+a^2)/2=a^2=3/20$.

In one dimension, the SAVE matrix is a scalar:

\widehat{\mathrm{SAVE}} =\frac12\left(1-\frac{27}{20}\right)^2 +\frac12\left(1-\frac{3}{20}\right)^2 =\frac12\frac{49}{400}+\frac12\frac{289}{400} =\frac{169}{400}=0.4225.

The example also shows why “deviation from identity” need not be a positive or negative scalar before squaring: the first slice variance exceeds one, while the second is below one. Both deviations contribute nonnegatively after squaring.

Aggregate the Squared Deviations

Once each $C_h$ has been calculated, form $I_p-C_h$, multiply that matrix by itself, and take the weighted sum. The square is essential: it records deviations in covariance without allowing positive and negative deviations from different slices to cancel as they might in a simple average.

Worked Example: A Diagonal SAVE Matrix and Its Leading Direction

Suppose three equally weighted slices have the following standardized covariance estimates: $C_1=\operatorname{diag}(0.2,0.3)$, $C_2=\operatorname{diag}(0.5,0.6)$, and $C_3=\operatorname{diag}(0.8,0.6)$. Calculate the SAVE matrix and its leading direction.

Because the matrices are diagonal, subtracting each from $I_2$ and squaring gives $(I_2-C_1)^2=\operatorname{diag}(0.64,0.49)$, $(I_2-C_2)^2=\operatorname{diag}(0.25,0.16)$, and $(I_2-C_3)^2=\operatorname{diag}(0.04,0.16)$. The weights are each $1/3$, so

\widehat{\mathrm{SAVE}} =\frac13 \begin{pmatrix} 0.64+0.25+0.04&0\\ 0&0.49+0.16+0.16 \end{pmatrix} = \begin{pmatrix} 0.31&0\\ 0&0.27 \end{pmatrix}.

The eigenvalues are $0.31$ and $0.27$, with corresponding unit eigenvectors $e_1=(1,0)^T$ and $e_2=(0,1)^T$. If one direction is retained, the leading estimated direction is $e_1$. This result comes from the larger weighted squared covariance deviation in the first coordinate.

Worked Example: Off-Diagonal Covariances and Eigenvector Extraction

Suppose two equally weighted slices have covariance estimates

C_1=\begin{pmatrix}0.2&0.2\\0.2&0.5\end{pmatrix}, \qquad C_2=\begin{pmatrix}1.3&-0.7\\-0.7&1.0\end{pmatrix}.

Both matrices are positive definite: their diagonal entries are positive, and their determinants are $0.2(0.5)-0.2^2=0.06$ and $1.3(1.0)-0.7^2=0.81$, respectively. The deviations from identity are

I_2-C_1=\begin{pmatrix}0.8&-0.2\\-0.2&0.5\end{pmatrix}, \qquad I_2-C_2=\begin{pmatrix}-0.3&0.7\\0.7&0\end{pmatrix}.

Multiplying each matrix by itself gives

(I_2-C_1)^2= \begin{pmatrix}0.68&-0.26\\-0.26&0.29\end{pmatrix}, \qquad (I_2-C_2)^2= \begin{pmatrix}0.58&-0.21\\-0.21&0.49\end{pmatrix}.

Therefore $\widehat{\mathrm{SAVE}}=\begin{pmatrix}0.63&-0.235\\-0.235&0.39\end{pmatrix}$. For a symmetric matrix $\begin{pmatrix}a&b\\b&d\end{pmatrix}$, the eigenvalues are $(a+d)/2\pm\sqrt{((a-d)/2)^2+b^2}$. Here the midpoint is $0.51$, and the square-root term is $\sqrt{0.12^2+0.235^2}=\sqrt{0.069625}\approx0.26387$. Thus the eigenvalues are approximately $0.77387$ and $0.24613$.

For the leading eigenvalue, the first row of the eigenvector equation gives $(0.63-0.77387)v_1-0.235v_2=0$, so $v_2\approx-0.6122v_1$. Normalizing gives a leading unit eigenvector approximately $(0.853,-0.522)^T$. Its sign could be reversed without changing the direction.

Why the Eigenvectors Give the Estimated Directions

The key algebraic property is that each squared deviation contributes a nonnegative amount in every direction. For any vector $v\in\mathbb{R}^p$,

v^T\widehat{\mathrm{SAVE}}v =\sum_{h=1}^H p_h\,v^T(I_p-C_h)^2v =\sum_{h=1}^H p_h\,\|(I_p-C_h)v\|^2\geq 0.

The second equality holds because $I_p-C_h$ is symmetric, so $v^T(I_p-C_h)^2v=\{(I_p-C_h)v\}^T\{(I_p-C_h)v\}$. This also proves that the SAVE estimate is positive semidefinite.

Theorem: The sample SAVE matrix is symmetric and positive semidefinite. For a chosen target dimension $d$, its leading $d$ orthonormal eigenvectors maximize the total quadratic form $ \sum_{j=1}^d v_j^T\widehat{\mathrm{SAVE}}v_j $ over all orthonormal $d$-tuples $v_1,\ldots,v_d$.

To verify the maximizing statement, let the eigenvalues be ordered $\lambda_1\geq\cdots\geq\lambda_p\geq0$, with orthonormal eigenvectors $u_1,\ldots,u_p$. Write any orthonormal $d$-tuple as $v_j=\sum_{k=1}^p a_{kj}u_k$, and let $q_k=\sum_{j=1}^d a_{kj}^2$. Orthonormality implies $0\leq q_k\leq1$ and $\sum_k q_k=d$. The total quadratic form is $\sum_k\lambda_kq_k$, which is at most $\sum_{k=1}^d\lambda_k$: weight assigned to an eigenvalue below $\lambda_d$ cannot improve on assigning that weight to an omitted eigenvalue at least $\lambda_d$. Taking $v_j=u_j$ for $j=1,\ldots,d$ attains this upper bound. This proves the claim. $\square$

These eigenvectors are expressed in standardized predictor coordinates. If $\hat{\eta}$ is a retained unit eigenvector, its direction on the original predictor scale is $\hat{\beta}=S^{-1/2}\hat{\eta}$, since $\hat{\eta}^Tz_i=\hat{\eta}^TS^{-1/2}(x_i-\bar{x})=\hat{\beta}^T(x_i-\bar{x})$. The resulting original-scale directions are orthonormal under the sample covariance metric: if the standardized eigenvectors form a matrix $U$, then $(S^{-1/2}U)^TS(S^{-1/2}U)=I$.

Standardize the predictors.
Compute $\bar{x}$ and $S$, then form $z_i=S^{-1/2}(x_i-\bar{x})$.

Estimate within-slice covariances.
For each response slice, calculate $\bar{z}_h$ and $C_h$ using the stated divisor.

Aggregate covariance deviations.
Compute $\sum_h(n_h/n)(I_p-C_h)^2$, using matrix multiplication for every square.

Extract and transform directions.
Order the eigenvectors by decreasing eigenvalue; transform retained directions by $S^{-1/2}$ if reporting them on the original scale.

Common Pitfalls and Practical Interpretation

A common mistake is to average the matrices $I_p-C_h$ first and then square the average. That computes $\{\sum_h p_h(I_p-C_h)\}^2$, which is generally not $\sum_h p_h(I_p-C_h)^2$. The first calculation can cancel deviations across slices; SAVE squares each slice's deviation before aggregation. Another error is squaring each matrix entry separately. The off-diagonal entries in the worked example show why the correct operation is matrix multiplication: the off-diagonal element of $A^2$ is $A_{11}A_{12}+A_{12}A_{22}$, not $A_{12}^2$.

Slice counts matter because a covariance estimate from a very small slice can be unstable, and a slice with only one observation has zero empirical covariance under the divisor-$n_h$ convention rather than a usable estimate of within-slice variation. Very unequal counts also make the weights unequal, as they should: the estimate targets the observed response distribution. Slice boundaries should be chosen with both response coverage and enough observations per slice in mind.

Key takeaway: SAVE estimation is a sequence of distinct operations: standardize globally, estimate covariance within each response slice, square each slice's deviation from identity, weight by slice proportion, and extract the leading eigenvectors. Keep the matrix order and the covariance divisor explicit.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Why does standardizing with $S^{-1/2}$ make the overall sample covariance equal to $I_p$?
If two slices have different counts, how are their weights determined, and what do those weights sum to?
What is the distinction between $\sum_h p_h(I_p-C_h)^2$ and $\{\sum_h p_h(I_p-C_h)\}^2$?
Why is the sample SAVE matrix positive semidefinite even when some entries of $I_p-C_h$ are negative?
If an eigenvector is estimated in standardized coordinates, what transformation expresses its direction on the original predictor scale?