<figure><img src="image_1.png" /></figure>

HAL
open science

# From explained variance of correlated components to PCA without orthogonality constraints

Marie Chavent, Guy Chavent

► To cite this version:

Marie Chavent, Guy Chavent. From explained variance of correlated components to PCA without orthogonality constraints. 2024. hal-04442489

HAL Id: hal-04442489

https://inria.hal.science/hal-04442489v1

Preprint submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# From explained variance of correlated components to PCA without orthogonality constraints

Marie Chavent *

Guy Chavent †

## Abstract

Block Principal Component Analysis (Block PCA) of a data matrix $A$, where loadings $Z$ are determined by maximization of $\|AZ\|^2$ over unit norm orthogonal loadings, is difficult to use for the design of sparse PCA by $\ell^1$ regularization, due to the difficulty of taking care of both the orthogonality constraint on loadings and the non differentiable $\ell^1$ penalty. Our objective in this paper is to relax the orthogonality constraint on loadings by introducing new objective functions $\text{expvar}(Y)$ which measure the part of the variance of the data matrix $A$ explained by correlated components $Y = AZ$. So we propose first a comprehensive study of mathematical and numerical properties of $\text{expvar}(Y)$ for two existing definitions Zou et al. [2006], Shen and Huang [2008] and four new definitions. Then we show that only two of these explained variance are fit to use as objective function in block PCA formulations for $A$ rid of orthogonality constraints.

**Keywords:** PCA, sparsity, dimension reduction, explained variance, orthogonality constraints, block optimization.

<footnote>* Université de Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, Talence, France
e-mail: marie.chavent@u-bordeaux.fr (corresponding author)</footnote>

<footnote>† Retired, collaborator to the SERENA Project Team, INRIA-Paris 2 rue Simone Iff, 75589 Paris, France</footnote>
# 1 Introduction

Many linear data analysis methods construct new variables that “best” summarize the columns of a $n \times p$ data matrix A where n observations are described by p numerical variables. The $m < p$ new variables are the columns of $Y = AZ$ where Z is the $p \times m$ loading matrix. These new variables are for instance principal components in PCA (Principal Component Analysis), canonical components in CA (Canonical Analysis) or PLS components in PLS (Partial Least Squares) regression. When the components are orthogonal (which is the case for PCA), it is usual to assess the quality of the components by measuring the part of the variance of A explained by $Y = AZ$ with :

$$ \expvar(AZ) = \|AZ\|_F^2, $$

where the subscript F denotes the Frobenius norm.

But this definition fails when the components are correlated. Correlated components appear when sparsity is introduced into the loading matrix Z to select the important original variables. The loading vectors and the components are no longer necessarily orthogonal. Two definitions have already been proposed to measure the explained variance of correlated components : the adjusted variance of Zou et al. [2006] and the total variance of Shen and Huang [2008]. Because there is no single definition, we introduce a set of three conditions to be satisfied by any explained variance definition and we propose a comprehensive study of mathematical and numerical properties of these two existing definitions together with four new definitions.

We prove first that the total variance of Shen and Huang [2008] called subspace explained variance hereafter, satisfies only two of three conditions. This lead us to propose other definitions by projection or by normalization on a set $X = [x_1, ..., x_m]$ of orthonormal vectors which “points in the same direction” as the components $Y = [y_1, ..., y_m]$. Such orthonormal vectors X can be associated to Y by three different rules: QR or UP (polar) decomposition of Y, or maximization of the explained variance. This leads to five definitions ( two normalized and three projected explained variance) :

- the QR and the UP normalized explained variances (QRnormVar and UPnormVar),

- the QR projected explained variance (QRprojVar) which is the adjusted variance of Zou et al. [2006],

- the UP projected explained variance (UPprojVar),

- the optimal projected explained variance (optprojVar).

We prove that the five above definitions satisfy also the two first compatibility conditions, but that only the three projected variances, satisfy the third one, and hence are proper explained variances.
Then we investigate both theoretically and numerically the size of the six explained variances, the existence of order between them, and their ability to rank in the same order - or not - different sets of components $Y$, which is for what they have been introduced at the origin.

Finally, we study the ability of the three proper definitions of explained variance to provide a substitute to the classical Block PCA formulation :

$$ \text{maximizing } \|AZ\|_F^2 \text{ under the orthonormality constraint } Z^T Z = I_m, \quad (1) $$

by replacing (1) by a *Explained Variance block PCA formulation* :

$$ \max_{\|z_j\|=1, j=1...m} \expvar(AZ), \quad (2) $$

which is rid of orthogonality constraints on loadings $Z = [z_1...z_m]$, and hence particularly suited as starting point for the design of sparse PCA algorithms.

As a conclusion of this study, we single out the *optimal projected explained variance* of $Y = AZ$:

$$ \expvar(Y) = \max_{X^T X = I_m} \sum_{j=1...m} \langle y_j, x_j \rangle^2, \quad (3) $$

which possesses, after introduction of weights, the desired properties in the sense that it is easy to compute and differentiate, and admits as unique maximizer the SVD solution of PCA $Z^* = [v_1...v_m]$ made of the $m$ first right singular vectors of $A$. This will be the starting point for the design of an efficient group-sparse block PCA algorithm in the companion paper [Chavent and Chavent, 2023].

The paper is organized as follows: In Section 2 we motivate the need for an explained variance definition, and establish necessary conditions to be satisfied by any definition of the variance explained by a set of non-necessarily orthogonal components. Section 3 establishes the properties of the subspace explained variance (total variance) of Shen and Huang [2008]. Section 4 is devoted to the definition of three projected (including the adjusted variance of Zou et al. [2006] and the optimal variance (3)) and two normalized explained variances, together with the study of their mathematical properties. Sections 5 and 5.3 present numerical comparisons of the size of the variance explained by the six definitions, and of their ability to rank in the same order (or not) non necessarily orthogonal components. Finally, Section 6 compares the ability of the three projected variance definitions to discriminate the singular value solution of PCA as their unique maximizer, which makes them fit to use as Block PCA objective function.

It should be noted that we have implemented the six explained variance definitions in R [R Core Team, 2021] in the package **sparsePCA** available at https://github.com/chavent/sparsePCA.
## 2 Defining variance explained by components when loadings are non orthogonal

We set us from now on in the context of Principal Component Analysis (PCA), where one seeks a small number $m$ of uncorrelated components $y_j = Az_j, j = 1...m$ by combining the $p$ columns of a data matrix $A$, each containing $n$ samples of a centered variable, with $m$ unit norm loading vectors $z_j, j = 1...m$, in such a way that the components $y_j$ retain the largest possible part of the total variance of $A$. We denote by

$$ \expvar(Y) = \text{the part of the variance of } A \text{ explained by } Y = AZ \quad (4) $$

and by

$$ \|A\|_F^2 = \sum_{i,j=1...m} a_{i,j}^2 \text{ the (total) variance of } A, \text{ where } \|\cdot\|_F \text{ is the Frobenius norm of } A. $$

At this point, $\expvar(Y)$ is a still loosely defined quantity outside of the solution of PCA given by :

$$ Z = V_m, \quad Y = U_m \Sigma_m = AV_m, \quad (5) $$

where $U_m$ and $V_m$ contain the $m$ first left and right singular vectors and $\Sigma_m$ is the diagonal matrix of the $m$ first singular values of the singular value decomposition (SVD) of $A$:

$$ A = U \Sigma V^T \quad \text{with} \quad U^T U = I_r, \quad V^T V = I_r, \\ \Sigma = \text{diag}(\sigma_1, \dots, \sigma_r) = r \times r \text{ matrix with } \sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0, $$

where $r$ is the rank of $A$, and the columns $u_1...u_r$ of $U$ and $v_1...v_r$ of $V$ are the left and right singular vectors of $A$.

The principal components $y_j = Az_j = \sigma_j u_j$ are orthogonal, and hence uncorrelated, so the sum of their variance $\|y_j\|^2$ represents the part of the total variance $\|A\|_F^2$ of $A$ explained by these $m$ principal components. So we see that the variance explained by the principal components $Y$ is:

$$ \expvar(Y) = \|Y\|_F^2 = \sum_{j=1...m} \|y_j\|^2 = \sum_{j=1...m} \sigma_j^2 \le \sum_{j=1...r} \sigma_j^2 = \|A\|_F^2. \quad (6) $$

However, when sparsity constraints are introduced into loading vectors like in sparse PCA for instance, non orthogonal loadings and components are generated and one has to face the problem of defining $\expvar(Y)$ for possibly non orthogonal components. Alas, formula (6) for $\expvar(Y)$ is strictly limited to the case of orthogonal components and loadings, as we see now.

Consider first the case of non orthogonal components: take for example for $Z$ an orthonormal basis of $\text{span}\{V_m\}$ but different from $V_m$. Then components $Y = AZ$ are not anymore orthogonal, and hence correlated, so the sum of their variances (the total variance of $Y$) is too optimistic, and one expects that

$$ \expvar(Y) < \|Y\|_F^2, \quad (7) $$
which shows that $||Y||_F^2$ is not a satisfying definition of $\expvar(Y)$ in this case.

Then take for example a data matrix $A$ with three singular values 3, 2, 1 and hence a total variance of 14. Then chose for $Z$ two linearly independent unit vectors close to the first right singular vector $v_1$. Then:

$$ ||AZ||_F^2 = ||Y||_F^2 = ||y_1||^2 + ||y_2||^2 = \underbrace{||Az_1||^2}_{\approx \sigma_1^2=9} + \underbrace{||Az_2||^2}_{\approx \sigma_1^2=9} \approx 18 > \underbrace{9+4+1}_{\sigma_1^2+\sigma_2^2+\sigma_3^2} = 14 = ||A||_F^2, $$

which violates property $\expvar(Y) \le ||A||_F^2$ implied by (4), so once again $||Y||_F^2$ is not suited as a definition of $\expvar(Y)$.

A last motivation for the search of definitions of $\expvar(Y)$: rather than solving the classical Block PCA formulation (1), why not solve the explained variance Block PCA formulation (2) by maximizing $\expvar(Y)$ under the sole unit norm constraint on the loadings $Z$? Sparse block PCA formulations based on such explained variance objective function eliminate the difficulty caused by the orthogonality constraints on the loadings, and, by construction, rule directly the balance between sparsity and explained variance.

So we address first in this paper the problem of defining the part of the variance of $A$ explained by components $Y = AZ$, under the sole condition that:

$$ ||z_j|| = 1, j = 1 \dots m, \quad \text{rank}(Z) = \text{rank}(Y) = m \le r. \quad (8) $$

In absence of a sound definition for the explained variance of correlated components, we define a set of hopefully reasonable necessary conditions to be satisfied by any such definition:

* **Condition 1**: when $Y, Z$ happen to be the SVD solution of PCA given by (5), $\expvar(Y)$ has to provide the *exact value* given by (6):

$$ \expvar(Y) = \sum_{j=1 \dots m} \sigma_j^2 \le \sum_{j=1 \dots r} \sigma_j^2 = ||A||_F^2 \quad \text{for } m = 1 \dots r. \quad (9) $$

* **Condition 2**: when $Y, Z$ satisfy only (8), the components are not orthogonal anymore, and one expects that, because of the correlation between the components, the part of the variance of $A$ explained by $Y$ is smaller than that of the PCA solution:

$$ \expvar(Y) \le \sum_{j=1 \dots m} \sigma_j^2. \quad (10) $$

* **Condition 3**: $\expvar(Y)$ has to take into account the possible correlation of the components as expected in (7):

$$ \expvar(Y) \le \sum_{j=1 \dots m} ||y_j||^2 = ||Y||_F^2, $$

with equality only if the components are orthogonal.
Any explained variance definition $expvar(AZ)$ which satisfies these conditions achieves its PCA maximum value $\sum_{j=1...r} \sigma_j^2$ (Condition 2) for $Z = V_m$ (Condition 1), and hence provides a block formulation for PCA without orthogonality constraints on loadings, which can be used to derive sparse PCA algorithms. This point of view will be developed in Section 6 below.

We propose now definitions for the explained variance of the components $Y = AZ$ associated to any matrix $Z$ of $m \le r$ linearly independent - but not necessarily orthogonal - unit norm loading vectors $z_j$, $j = 1, ..., m$. These definitions will include those introduced by Shen and Huang [2008] (see section 3 below) and Zou et al. [2006] (see section 4.2).

## 3 Subspace explained variance

We start here from a reformulation of the explained variance (6) of the principal components $Y = AZ$, based on the subspace spanned by $Z = V_m$. Let $P_{V_m}$ denotes the orthogonal projection on this subspace. Then $P_Z = P_{V_m} = V_m V_m^T$, so that:

$$ \expvar(Y) = \|Y\|_F^2 = \|AP_Z\|_F^2 = \|AP_{V_m}\|_F^2 = \sum_{j=1...m} \sigma_j^2. \quad (11) $$

When $Z$ satisfies (8) only, we proceed by analogy with (11), and define the subspace explained variance of $Y = AZ$ by:

$$ \expvar_{subsp}(Y) \stackrel{\text{def}}{=} \|AP_Z\|_F^2 = \text{tr}\{Y^T Y(Z^T Z)^{-1}\}, \quad (12) $$

where $P_Z = Z(Z^T Z)^{-1}Z^T$ denotes the projection matrix on the subspace spanned by $Z$. This shows that subspace explained variance is the Rayleigh quotient associated to $A, Z$. Note that with this definition, where the explained variance depends only of the subspace spanned by $[z_1...z_m]$, the normalization of the loadings $z_j$ is not necessary. Of course, we will still continue to represent loadings by unit norm vectors - but this is here only a convenience.

The subspace explained variance coincides with the total explained variance introduced by Shen and Huang [2008, section 2.3 p. 1021], which they proved was increasing with the number of loadings, and bounded by the variance $\|A\|_F^2$ of the data. The next lemma gives a complete picture of its properties:

**Lemma 3.1.** (Subspace explained Variance) Let $Z$ satisfy (8). Then the subspace explained variance defined by (12) satisfies:

$$ \expvar_{subsp}(Y) = \sigma_1^2 + \cdots + \sigma_m^2 = \text{Max} \quad \Leftrightarrow \quad \text{span}\{Z\} = \text{span}\{V_m\}. \quad (13) $$

and hence satisfies conditions 1 and 2. However, it does not satisfy condition 3:

- when the components $Y = AZ$ happen to be orthogonal:

$$ \expvar_{subsp}(Y) \ge \|Y\|_F^2 \text{ with equality iff } Z \stackrel{\text{perm}}{=} V_m, \quad (14) $$
where “perm” denotes the equality of matrices up to a column permutation and multiplicatipn by ±1,

- when the loadings Z happen to be orthogonal :

$$
\expvar_{subsp}(Y) = \|Y\|_F^2 \quad \text{with non necessarily orthogonal components.} \tag{15}
$$

The proof is in Section 8.2 of the Appendix. This lemma shows that $\expvar_{subsp}(Y)$
verifies only conditions 1 and 2, and overestimates the explained variance in two
cases :

- when components $Y$ are orthogonal without pointing in the direction of
the left singular vectors, inequality (14) implies $\expvar_{subsp}(Y) > \|Y\|_F^2$,
which contradicts condition 3,

- when loadings Z are orthogonal without pointing in the direction of the
right singular vectors, then the components are correlated and (15) contra-
dicts condition 3, which requires in that case that $\expvar_{subsp}(Y) < \|Y\|_F^2$.

We explore in the next section other directions in the hope of being able to
comply with all conditions 1, 2 and 3.

4 Projected and normalized explained variances

We start now from definition (6) of the explained variance in the case of PCA.
A natural generalization would be :

$$
\expvar(Y) \stackrel{?}{=} \sum_{j=1\dots m} \|y_j\|^2 = \|Y\|_F^2 = \|AZ\|_F^2,
$$

which, as we have seen in the Section 2, is not anymore an acceptable defi-
nition when components are correlated.

However, this definition continues to make perfect sense for the variance
explained by components as long as they are orthogonal, without pointing nec-
essarily in the direction of left singular vectors. Hence a natural way to eliminate
the redundancy caused by the orthogonality default of the components *Y* and
*to satisfy Condition 3 is to :*

1. **Step 1:** choose a rule to associate to the components *Y* an *orthonormal basis* *X* of span{*Y*} that, loosely speaking, “points in the direction of the components *Y*”, and, when the components *Y* happen to be orthogonal, points in the directions of *Y* itself. So the rule for the choice of the basis *X* associated to *Y* has to satisfy :

$$
\begin{align}
& X^T X = I_m, \quad \text{span}\{X\} = \text{span}\{Y\}, \tag{16} \\
& < y_j, x_j > \ge 0 \quad \forall j = 1, \dots, m, \nonumber \\
& \langle y_j, y_k \rangle = 0 \quad \forall j \ne k \quad \Longrightarrow \quad x_j = y_j / \|y_j\| \quad \forall j = 1, \dots, m. \nonumber
\end{align}
$$
Examples of such rules are :

QR decomposition : $Y \rightsquigarrow X = Q$ solution of $Y = QR, Q^TQ = I_m$ where $R$ is an upper triangular matrix,

Polar decomposition : $Y \rightsquigarrow X = U$ solution of $Y = UP, U^TU = I_m$ where $P^T = P \in \mathbb{R}^{m \times m}, P \ge 0$. (17)

2. **Step 2**: associate to $Y$ orthogonal adjusted components $Y'$ along the $X$ axes, and define the variance explained by the components $Y$ by :

$$ \expvar(Y) \stackrel{\text{def}}{=} \|Y'\|_F^2. $$

Two ways for obtaining the adjusted components $Y'$ are considered hereafter: by projection or by normalization as illustrated in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of projected and normalized explained variances : The ellipse $\mathcal{E}^2$ represents the image by $A$ of all unit norm loadings. Let the $X = [x_1, x_2]$ be the orthonormal vectors associated to the correlated components $Y = [y_1, y_2]$. Then the extremities of the adjusted components $Y' = [y'_1, y'_2]$ obtained by projection are the two blue dots on the $x_1, x_2$ axes, and those obtained by normalization are the two green dots located at the intersection of the axes with $\mathcal{E}^2$.</figcaption></figure>

## 4.1 Normalized explained variances

We consider first in this section the case where the adjusted components $Y'$ in step 2 are obtained by “normalization” in the directions of the chosen orthonormal basis $X$ of span{$Y$}. More precisely, the idea is to choose the abscissa of
$y'_j$ on the $x_j$ axis by requiring that $y'_j$ is the image by $A$ of some unit norm adjusted loading $z'_j$. This is illustrated in Figure 1: one associates to components $Y = [y_1, y_2]$ the adjusted components $Y' = [y'_1, y'_2]$ whose extremities are the points on the $x_1, x_2$ axes located on the ellipse image by $A$ of the unit ball of the loading space. In order to determine the adjusted loadings $z'_j$, one computes first (non necessarily unit norm) loadings $T = [t_1, \dots, t_m]$ by performing on the loadings $z_j$ the same linear combinations $M$:

$$Y = XM, \qquad (18)$$

that transformed $Y$ into $X$, which leads to define $T$ by:

$$Z = TM. \qquad (19)$$

Multiplying then by $A$ left gives $Y = AZ = ATM$ and, comparison with (18) shows that $X = AT$ and so $x_j = At_j, j = 1 \dots m$. The adjusted components $y'_j$ are then defined by:

$$y'_j = Az'_j \text{ with } z'_j = t_j/\|t_j\| \text{ so that } \|z'_j\| = 1, \quad y'_j = x_j/\|x_j\|, \quad j = 1 \dots m, \qquad (20)$$

and the normalized explained variance of $Y$ estimated with $X$ is defined by:

$$\expvar_{norm}^{X}(Y) = \|Y'\|_{F}^{2} = \sum_{j=1 \dots m} 1/\|t_j\|^{2}, \qquad (21)$$

where $T = [t_1, \dots, t_m]$ is given by (19) and (18). Before specifying the rule of step 1 which associates an orthonormal basis $X$ to the components $Y$, we give some properties of $\expvar_{norm}^{X}(Y)$ which hold independently of the chosen basis $X$.

**Lemma 4.1.** (Normalized explained variance) For any unit norm loadings $Z$ and any basis $X$ chosen according to the rule (16), the normalized explained variance of $Y = AZ$ defined by (21) satisfies conditions 1 and 2 and:

$$\expvar_{norm}^{X}(Y) \le \expvar_{subsp}(Y) \le \sigma_{1}^{2} + \cdots + \sigma_{m}^{2}. \qquad (22)$$

Lemma 4.1 follows from Lemma 3.1 applied to the orthogonal components $Y'$:

$$\expvar_{norm}^{X}(Y) \stackrel{\text{def}}{=} \|Y'\|_{F}^{2} \le \expvar_{subsp}(Y') \le \sigma_{1}^{2} + \cdots + \sigma_{m}^{2}, \qquad (23)$$

which proves (22). We can now specify the rules for the choice of $X$ to define two normalized explained variance satisfying conditions 1 and 2.

**QR normalized variance.** Let $X$ be defined by the QR-decomposition $Y = XR$ of the components $Y = AZ$ as recalled in (17). Then (21) leads to another definition of variance:

$$\expvar_{norm}^{QR}(Y) = \sum_{j=1 \dots m} 1/\|t_j\|^2, \quad Z = TR, \ R \text{ upper triangular.} \qquad (24)$$

The normalized variance $\expvar_{norm}^{QR}(AZ)$ does not satisfy condition 3, as the counter example of Figure 2, left, shows.
<figure><img src="image_4.png" /><figcaption>Figure 2: Counterexamples for property 3. **Left**: Let $X$ denote the basis associated to components $Y$ by the QR-decomposition, and $Y'$ be the corresponding normalized adjusted components. One sees that $\expvar_{norm}^{QR}(Y) = \|Y'\|^2 \ge \|Y\|^2$, which violates property 3. **Right**: Let $X$ denote the basis associated to components $Y$ by the polar decomposition, and $Y'$ be the corresponding normalized adjusted components. One sees that $\expvar_{norm}^{QR}(Y) = \|Y'\|^2 \ge \|Y\|^2$, which violates property 3.</figcaption></figure>

**UP normalized variance.** Let $X$ be defined by the UP-decomposition (polar decomposition) $Y = XP$ of the components $Y = AZ$ as recalled in (17). Then (21) defines another variance :

$$ \expvar_{norm}^{UP}(Y) = \sum_{j=1...m} 1/\|t_j\|^2, \quad Z = TP, \quad P = (Y^T Y)^{1/2}. \quad (25) $$

The normalized variance $\expvar_{norm}^{UP}(AZ)$ does not either satisfy condition 3, as the counter example of Figure 2, right, shows.

So we explore in the next section another road in order to comply with all conditions 1,2 and 3.

## 4.2 Projected explained variances

We consider in this section the case where the adjusted components $Y'$ in step 2 are obtained by projection of the components $Y$ on the chosen orthonormal basis $X$ of span{$Y$}. The adjusted component $y'_j$ is hence defined by :

$$ y'_{j} = \langle y_{j}, x_{j} \rangle x_{j}, \quad j = 1 \dots m. $$

The so called projected explained variance of $Y$ estimated with $X$ is then defined by :

$$ \expvar_{proj}^{X}(Y) \stackrel{\text{def}}{=} \|Y'\|_{F}^{2} = \sum_{j=1...m} \langle y_{j}, x_{j} \rangle^{2} \quad (26) $$
Before specifying the rule of step 1 which associates an orthonormal basis X to the possibly correlated components Y, we give some properties of $expvar_{proj}^X(Y)$ which hold independently of the chosen basis X:

**Lemma 4.2.** (Projected explained variances) For any unit norm loadings Z and any basis X chosen according to the rule (16), the projected explained variance of Y = AZ defined by (26) satisfies conditions 1, 2 and 3 and:

$$ \expvar_{proj}^X(Y) \leq \expvar_{subsp}(Y) \leq \sigma_1^2 + \dots + \sigma_m^2. \quad (27) $$

The proof is in Section 8.3 of the Appendix. We can now specify the rules for the selection of the orthonormal basis X, which give each, according to Lemma 4.2, a projected explained variance satisfying conditions 1 to 3.

**QR projected explained variance.** We choose here to associate to the components Y, the orthonormal basis X = Q obtained by QR-decomposition of Y as recalled in (17). The vector $x_1$ of the basis X = Q is chosen in the direction of the component with the larger norm, and the remaining components are projected on the orthogonal subspace to $x_1$. Then $x_2$ is determined by the same process applied in the orthogonal subspace, and so on. This reordering ensures that the basis X = Q associated to Y will point in the direction of Y primarily for the components of larger variance. The (order dependent) resulting QR projected explained variance is given by:

$$ \expvar_{proj}^{QR}(Y) = \sum_{j=1 \dots m} \langle y_j, x_j \rangle^2 = \sum_{j=1 \dots m} r_{j,j}^2. \quad (28) $$

It coincides with the adjusted variance introduced in Zou et al. [2006].

**UP Projected explained variance.** We choose now to associate to the components Y, the orthonormal basis X = U obtained by UP-decomposition (polar decomposition) of Y as recalled in (17). The basis X = U does its best to point in the same direction as the components Y, in that it maximizes the scalar product $\langle Y, X \rangle_F = \sum_{j=1 \dots m} \langle x_j, y_j \rangle$. The (order independent) resulting UP projected explained variance is given by:

$$ \expvar_{proj}^{UP}(Y) = \sum_{j=1 \dots m} \langle y_j, x_j \rangle^2 = \sum_{j=1 \dots m} p_{j,j}^2. \quad (29) $$

**Optimal projected explained variance.** The idea here is to associate to the components Y, the basis X which gives the largest projected explained variance defined in (26). This choice satisfies obviously condition (16) and the so called optimal projected explained variance is defined by:

$$ \expvar_{proj}^{opt}(Y) = \max_{X^T X = I_m} \sum_{j=1 \dots m} \langle y_j, x_j \rangle^2, \quad (30) $$
The numerical computation of the optimal projected explained variance requires the maximization of the convex function $X \rightsquigarrow \sum_{j=1...m} \langle y_j, x_j \rangle^2$ under the constraint $X^T X = I_m$. This can be done using the algorithm of Journée et al. [2010], which gives here :

$$X_{k+1} = \operatorname{polar}(2Y \operatorname{diag}(X_k^T Y)), \quad X_0 = U = \operatorname{polar}(Y). \qquad (31)$$

where $2Y \operatorname{diag}(X_k^T Y)$ is the gradient at $X_k$ of the function $X \rightsquigarrow \sum_{j=1...m} \langle y_j, x_j \rangle^2$.

# 5 Comparison of the explained variances

In this section, the six variances explained by non orthogonal components (see Table 1), are compared theoretically and numerically.

<table><caption>Table 1: Summary of the 6 variance definitions.</caption><thead><tr><th>Name</th><th>Notation</th><th>Definition</th><th>Short name</th></tr></thead><tbody><tr><td>subspace variance</td><td>expvar<sub>subsp</sub></td><td>[12]</td><td>subspVar</td></tr><tr><td>QR normalized variance</td><td>expvar<sup>QR</sup><sub>norm</sub></td><td>[24]</td><td>QRnormVar</td></tr><tr><td>UP normalized variance</td><td>expvar<sup>UP</sup><sub>norm</sub></td><td>[25]</td><td>UPnormVar</td></tr><tr><td>QR projected variance</td><td>expvar<sup>QR</sup><sub>proj</sub></td><td>[28]</td><td>QRprojVar</td></tr><tr><td>UP projected variance</td><td>expvar<sup>UP</sup><sub>proj</sub></td><td>[29]</td><td>UPprojVar</td></tr><tr><td>optimal projected variance</td><td>expvar<sup>opt</sup><sub>proj</sub></td><td>[30]</td><td>optprojVar</td></tr></tbody></table>

## 5.1 What we know

We give first theoretical results on the relative magnitudes of the 6 explained variances :

- The subspace explained variance is larger than any of the five other variances (Lemmas 4.2 and 4.1). It is even larger than expected when components or loadings are orthogonal without being left and right singular vectors (see (14) and (15)).

- The optimal projected variance is greater by definition than any other projected variance, in particular greater than the QR and the UP projected variance.

- There is no natural order between the QR and UP projected variances : when the components $Y$ are of equal norm, the basis $X$ which maximizes $\expvar_{proj}^X(Y)$ is $\operatorname{polar}(Y)$, which implies in particular that :

$$\expvar_{proj}^{opt}(Y) = \expvar_{proj}^{UP}(Y) \geq \expvar_{proj}^{QR}(Y). \qquad (32)$$

But the converse of the last inequality can hold when the norms of the components are very different : for $m = 2$, one checks that $\|y_2\|/\|y_1\|$ small enough implies that $\expvar_{proj}^{UP}(Y) \leq \expvar_{proj}^{QR}(Y)$.
- There is no natural order between the QR and UP normalized variances : for components Y such that the basis X associated by QR-decomposition coincides with the m-first left singular vectors $U_m$ of A, one has, according to Lemma 4.1 :

$$ \expvar_{norm}^{QR}(Y) = \sum_{j=1...m} \sigma_j^2 \geq \expvar_{norm}^{UP}(Y), $$

with a strict inequality as soon as Y and $U_m \text{diag}\{\sigma_j\}$ don't coincide. The same reasoning with the polar decomposition in place of the QR decomposition shows that the converse inequality can happen.

- There is no natural order between the variances defined by projection and normalization, as illustrated in Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: The two sets of components $Y = [y_1, y_2]$ and $\tilde{Y} = [\tilde{y}_1, \tilde{y}_2]$ have been chosen such that their polar decomposition produces the same basis $X = [x_1, x_2]$, and one sees that : $\expvar_{proj}^{UP}(\tilde{Y}) \leq \expvar_{norm}^{UP}(\tilde{Y}) = \expvar_{norm}^{UP}(Y) \leq \expvar_{proj}^{UP}(Y)$.</figcaption></figure>

## 5.2 What we see

We compare now numerically the six explained variances. The comparison is made on non orthogonal components Y obtained by applying sparse PCA to simulated data matrices A. The matrices A are obtained using the simulation scheme of Chavent and Chavent [2023] based on m = 4 underlying loadings vectors of size p = 20 and m = 4 underlying first eigenvalues which are chosen to be either close or different. More precisely, 100 matrices A of size n × 20 were drawn randomly using either the “close eigenvalues” or the “different eigenvalues” scheme. Three sparse PCA algorithms [see Chavent and Chavent, 2023]
were applied to each matrix A for a grid of 101 values of sparsity parameters $\lambda \in [0, 1]$. Finally, $30300 = 100 * 3 * 101$ loadings matrices Z and components matrices Y = AZ where obtained for the “close eigenvalues” and for the “different eigenvalues” scheme. The variance explained by these components where performed using the six variance definitions to finally obtain the dimensionless proportion of explained variance (pev) defined by:

$$pev = \frac{\expvar(Y)}{\|A\|_F^2} \le \frac{\sigma_1^2 + \cdots + \sigma_m^2}{\|A\|_F^2},$$

where the right inequality follows from (10) in Condition 2, which is satisfied by all definitions.

Figure 4 gives the mean pev for the 300 non orthogonal components (for the “close eigenvalues” and the “different eigenvalues” case) as a function of the sparsity parameter $\lambda$ and for each definition of explained variance.

<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison of the mean pev (proportion of explained variance) over two sets of components (top and bottom) as function of sparsity parameter $\lambda$ for the six variance definitions.</figcaption></figure>

For the “close eigenvalues” case (top), the six definitions give relatively close results and the results produced by optprojVar and UPprojVar are so close that
they cannot be distinguished on the figure (remember that they would coincide were the norms equal, see (32)). One sees also that the results seem to be in a certain order for all $\lambda$:

$$ \text{subspVar} \geq \text{optprojVar} \geq \text{UPprojVar} \geq \text{QRprojVar} \geq \text{QRnormVar} \geq \text{UPnormVar} $$

(the two first inequalities are not a surprise, as they hold theoretically). For the “different eigenvalues” case (bottom), a zoom on the curves shows that subspVar and optprojVar are again larger than all other variances. But the previous apparent order between UPprojVar, QRprojVar and QRnormVar is not longer observed. The UPnormVar remains the smallest but its behavior seems disturbed as the sparsity parameter $\lambda$ increases.

Table 2 shows that all definitions but UPnormVar (in the “different eigenvalues” case) exhibit quite similar dispersions over the 300 realizations (for $\lambda = 0.3$).

<table><caption>Table 2: Standard deviations×100 of the six pev (proportion of explained variance) obtained for $\lambda = 0.3$ with the three algorithms over the two sets of components (close eigenvalues and different eigenvalues).</caption><thead><tr><th></th><th>Close eigenvalues</th><th>Different eigenvalues</th></tr></thead><tbody><tr><td>subspVar</td><td>0.63</td><td>1.63</td></tr><tr><td>optprojVar</td><td>0.66</td><td>1.38</td></tr><tr><td>UPprojVar</td><td>0.66</td><td>1.06</td></tr><tr><td>QRprojVar</td><td>0.72</td><td>1.25</td></tr><tr><td>QRnormVar</td><td>1.06</td><td>1.33</td></tr><tr><td>UPnormVar</td><td>1.21</td><td>7.74</td></tr></tbody></table>

## 5.3 Ranking properties

The proportions of explained variance are meant to be used for the ranking of algorithms, so it is important to figure out whether or not definitions $i$ and $j$ of explained variance will rank in the same order the components $Y$ and $Y'$ obtained from possibly different algorithms and/or sparsity parameter $\lambda$ and/or realization of the data matrix $A$. The components obtained with the 3 algorithms, the 50 smallest values $\lambda$ and 100 realizations of $A$ (“different eigenvalues” case) gave 15000 × 14999/2 couples $(Y, Y')$ to be tested. Among these couples, we may consider as $\epsilon$-distinguishable from the point of view of our explained variances those for which

$$ |\text{pev}_i(Y) - \text{pev}_i(Y')| \geq \epsilon \quad \text{for all} \quad i = 1 \dots 6 $$

for some $\epsilon \geq 0$. Table 3 shows the percentage of cases where $\text{pev}_i$ and $\text{pev}_j$ rank identically components $Y$ and $Y'$ among all $\epsilon$-distinguishable couples.

The good news is that all three projected variances $\text{optprojVar}$, $\text{UPprojVar}$ and $\text{QRprojVar}$, as well as the normalized variance $\text{QRnormVar}$, produce the
<table><caption>Table 3: The entry of each table on line *i* and column *j* gives the percentage of $\epsilon$-distinguishable couples $Y, Y'$ which are ranked identically by $\text{pev}_i$ and $\text{pev}_j$. Top: $\epsilon = 0$, middle: $\epsilon = 10^{-3}$, bottom: $\epsilon = 10^{-2}$.</caption><thead><tr><th></th><th>optprojVar</th><th>UPprojVar</th><th>QRprojVar</th><th>QRnormVar</th><th>UPnormVar</th></tr></thead><tbody><tr><td>subspVar</td><td>79.98</td><td>71.05</td><td>70.71</td><td>69.75</td><td>56.19</td></tr><tr><td>optprojVar</td><td></td><td>88.93</td><td>89.87</td><td>89.19</td><td>73.71</td></tr><tr><td>UPprojVar</td><td></td><td></td><td>96.22</td><td>95.27</td><td>84.62</td></tr><tr><td>QRprojVar</td><td></td><td></td><td></td><td>98.34</td><td>83.15</td></tr><tr><td>QRnormVar</td><td></td><td></td><td></td><td></td><td>82.75</td></tr><tr><td>UPnormVar</td><td></td><td></td><td></td><td></td><td></td></tr></tbody><tbody><tr><td>subspVar</td><td>86.13</td><td>83.78</td><td>84.70</td><td>84.57</td><td>68.68</td></tr><tr><td>optprojVar</td><td></td><td>96.40</td><td>98.57</td><td>98.32</td><td>81.14</td></tr><tr><td>UPprojVar</td><td></td><td></td><td>97.81</td><td>97.72</td><td>84.74</td></tr><tr><td>QRprojVar</td><td></td><td></td><td></td><td>99.66</td><td>82.55</td></tr><tr><td>QRnormVar</td><td></td><td></td><td></td><td></td><td>82.66</td></tr><tr><td>UPnormVar</td><td></td><td></td><td></td><td></td><td></td></tr></tbody><tbody><tr><td>subspVar</td><td>89.57</td><td>89.57</td><td>89.57</td><td>89.57</td><td>68.80</td></tr><tr><td>optprojVar</td><td></td><td>100.00</td><td>100.00</td><td>100.00</td><td>79.23</td></tr><tr><td>UPprojVar</td><td></td><td></td><td>100.00</td><td>100.00</td><td>79.23</td></tr><tr><td>QRprojVar</td><td></td><td></td><td></td><td>100.00</td><td>79.23</td></tr><tr><td>QRnormVar</td><td></td><td></td><td></td><td></td><td>79.23</td></tr><tr><td>UPnormVar</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

same ranking as soon as one considers that differences in proportion of explained variance under $10^{-2}$ are not significative. For the same $\epsilon$, the two other definitions subspVar and UPnormVar still produce quite different rankings. The numerical results of this section show that all investigated definitions rank the explained variance of components in essentially the same order, and hence can all be used to compare the variance explained by components obtained by different algorithms or with different parameters. Of course, this is only an experimental result based on our simulated data sets, that needs to be confirmed by further numerical tests.

# 6 Explained variance block PCA formulations

With the objective in mind to get rid of the orthogonality constraints in the usual block PCA formulation (1), we discuss in this section the possibility of using the above defined explained variance measures as objective function in new block PCA formulations where the loadings are subject to unit norm constraints
only :

$$ \max_{\|z_j\| = 1, j = 1 \dots m} \expvar(AZ) = \sigma_1^2 + \dots + \sigma_m^2, \quad (33) $$

As seen in Lemmas 3.1, 4.1 and 4.2, property (33) holds for the subspace variance and the five normalized and projected variance definitions. However, from these six definitions, only the three projected explained variances satisfy (Lemma 4.2) all necessary conditions 1, 2 and 3 of section 2 required for a proper definition of explained variance. So we shall limit our search for an explained variance block formulation to these three definitions.

A projected explained variance $\expvar_{proj}(AZ)$ will provide a block PCA formulation (33) if and only if its sole maximizers are the loadings Z made of the m first right singular vectors $V_m$ in any order, up to a ±1 multiplication, or in short:

$$ Z \stackrel{\text{perm}}{=} V_m \text{ are the sole maximizers of } \expvar(AZ). \quad (34) $$

We denote also by

$$ A_m \stackrel{\text{def}}{=} U_m \Sigma_m V_m^T \quad \text{with} \quad \Sigma_m \stackrel{\text{def}}{=} \text{diag}\{\sigma_1 \dots \sigma_m\} \quad (35) $$

the restriction of the matrix A to the subspace of the m first right singular vectors $V_1 \dots V_m$. It is a bijection from the right to the left singular subspaces associated to the $m \le r$ first singular values. In the singular basis $U_m, V_m$, $A_m$ reduces to $\Sigma_m$.

The next lemma characterizes the maximizers of any projected variance:

**Lemma 6.1.** *Maximizers of projected variances*

1. For a given set of components Y, the basis X which gives the largest projected explained variance $\expvar_{proj}^X(Y)$ is necessarily a solution of:

$$ Y \text{diag}(X^T Y) = X P \text{ with } P^T = P \text{ and } P \ge 0 \quad (36) $$

2. The projected explained variance $\expvar_{proj}^X(AZ)$ achieves it PCA maximum value $\sigma_1^2 + \dots + \sigma_m^2$ if and only if Z and X satisfy the three following conditions:

$$ \text{span}\{Z\} = \text{span}\{V_m\} \quad (37) $$

the loadings Z are $(A_m^T A_m)^{-1}$ orthogonal (38)

the basis X associated to $Y = AZ$ is $X = [n_1 \dots n_m]$. (39)

where $n_1 \dots n_m$ denote unit normals at $y_1 \dots y_m$ to the ellipsoid $\mathcal{E}^m$ of $\text{span}\{U_m\}$ image by $A_m$ of the unit sphere of $\text{span}\{V_m\}$ (see Figure 5)
<figure><img src="image_7.png" /><figcaption>Figure 5: Illustration of condition (39) for $m = 2$. The ellipse $\mathcal{E}^2$ is the image by $A$ of all unit norm loadings $z$ of the right singular space $V_2$, with half axes $\sigma_1$ and $\sigma_2$. Let $y_1, y_2$ be two (possibly correlated) components such that the normals $n_1, n_2$ to $\mathcal{E}^2$ at $y_1, y_2$ are orthogonal, the tangents are orthogonal too and hence meet on the Cartan Circle of radius $R = (\sigma_1^2 + \sigma_2^2)^{1/2}$. Then for the choice $X = (n_1, n_2)$ one has $\expvar_{proj}^X Y = \|y'_1\|^2 + \|y'_2\|^2 = R^2 = \sigma_1^2 + \sigma_2^2$.</figcaption></figure>

The proof is in Appendix 8.4

Condition (38) reduces, in the singular basis $V_m$, to $\Sigma_m^{-2}$-orthogonality of the loadings. It ensures that the normals $n_j$ are orthogonal, and hence can be chosen as the basis $X$ associated to $Y$ in (39). Of course, the SVD solutions $Z \stackrel{\text{perm}}{=} V_m$ of PCA satisfy always (37) (38) (39).

Once a rule $Y \rightsquigarrow X$ has been specified, points 1 and 2 of the lemma will make it possible to determine whether or not (34) is satisfied.

## 6.1 Maximizers of $\expvar_{proj}^{QR}(Y)$

The projected explained variance $\expvar_{proj}^{QR}(AZ)$ attains its maximum $\sigma_1^2 + \dots + \sigma_m^2$ if and only if:

$$Z \stackrel{\text{perm}}{=} V_m, \qquad (40)$$

When the maximum norm selection procedure is applied at each step and the components renumbered accordingly, the unique maximizer is $Z = V_m$.

**Proof:** The “if” part of (40) is trivial, we prove the “only if” part: let $Y$ be such that $\expvar_{proj}^{QR} Y = \sigma_1^2 + \dots + \sigma_m^2$, and $X$ be given by its QR decomposition
$Y = XR$ with $R$ upper triangular. By hypothesis, $X$ maximizes $\expvar^X Y$, and point 1 of Lemmas 6.1 implies that $Y\mathrm{diag}(X^T Y) = XP = XR\mathrm{diag}(X^T Y)$, with $P = R\mathrm{diag}(X^T Y)$ symmetric positive. This implies that $R$ is a diagonal matrix, so that all $x_j$ point in the direction of $y_j$. But property 2 of Lemma 6.1 shows that the $x_j$'s are also normal to the ellipsoid $\mathcal{E}^m$ at $y_j$, which can happen only if $y_j = Az_j$ coincides with its principal axes and hence each $z_j$ is one of the $m$ first right eigenvectors $v_j$.

**Proposition 6.2.** The QR projected explained variance $\expvar_{proj}^{QR}(AZ)$ (the adjusted variance of Zou et al. [2006]) associated to unit norm loadings Z provides a block PCA formulation:

$$ \max_{\|z_j\| = 1, j = 1 \dots m} \expvar_{proj}^{QR}(AZ) = \sigma_1^2 + \dots + \sigma_m^2, \quad (41) $$

which admits the SVD solution (40) as “unique” maximizer.

Numerical implementation of this formulation requires the computation of the gradient of the $Z \rightsquigarrow \expvar_{proj}^{QR}(AZ)$ function, which is defined through the QR decomposition of $Y = AZ$. This can be done by the adjoint state method Chavent [2010], which is feasible but may be cumbersome. The block PCA formulation (41) can be used as starting point for the design of sparse PCA algorithms, keeping in mind that enforcing sparsity by subtracting the $\ell^1$ norm of loadings leads to a difficult, though tractable, non smooth optimization problem.

## 6.2 Maximizers of $\expvar_{proj}^{UP}(Y)$

The projected explained variance $\expvar_{proj}^{UP}(AZ)$ attains its maximum $\sigma_1^2 + \dots + \sigma_m^2$ if and only if:

$$ Z \stackrel{\text{perm}}{=} V_m \quad \text{(SVD solution) or} \quad Z = Z^\# \quad (42) $$

where the “parasitic” solution $Z^\#$ is such that the components $Y^\# = AZ^\#$ satisfy $\langle y_1^\#, x_1 \rangle = \dots = \langle y_m^\#, x_m \rangle$, with the hyperplanes tangent to $\mathcal{E}^m$ at $y_j^\#$ delimiting an $m$-dimensional orthant of $\mathrm{span}\{V_m\}$. For $m = 2$, the parasitic solution is illustrated on figure 6, where one sees that the components $y_j^\#$ (in red) corresponding to the choice of different principal axes for the intersection of the tangents coincide up to a multiplication by $\pm 1$; the SVD components are in blue. The proof is in Appendix 8.5.

So $\expvar_{proj}^{UP}(AZ)$ cannot be used for the construction of a block PCA formulation like (33), as the optimization algorithm might converge to the parasitic solution $Z^\#$!
<figure><img src="image_8.png" /><figcaption>Figure 6: Illustration, for $m = 2$, of the parasitic maximizer $Y^\#$ of $\text{expvar}_{\text{proj}}^{\text{UP}}$. The polar decomposition associates to components $Y^\# = y_1^\#, y_2^\#$ orthonormal vectors $X = x_1, x_2$ such that $\|y_1^\#\| = \|y_2^\#\|$. The red points on the Cartan circle correspond to the parasitic maximizers $Y^\# = \pm y_1^\#, \pm y_2^\#$, the blue points to the SVD solution $Y^* = \pm y_1^*, \pm y_2^* = \pm \sigma_1 u_1, \pm \sigma_2 u_2$</figcaption></figure>

## 6.3 Maximizers of $\text{expvar}_{\text{proj}}^{\text{opt}}(Y)$

According to point 2 of Lemma 6.1, the optimal projected explained variance $\text{expvar}_{\text{proj}}^{\text{opt}}(AZ)$ attains its maximum $\sigma_1^2 + \cdots + \sigma_m^2$ if and only if :

$$ \text{span}\{Z\} = \text{span}\{V_m\} \quad \text{and} \quad \text{the loadings } Z \text{ are } (A_m^T A_m)^{-1} \text{ orthogonal.} \tag{43} $$

This situation is similar to the maximization (1) of $\|AZ\|^2$ under orthonormality constraints $Z^TZ = I_m$, were the maximum is attained for all orthonormal $Z$ which $\text{span} V_m$. But the difference is that the orthogonality condition (43)-right is not a constraint for the maximization of $\text{expvar}_{\text{proj}}^{\text{opt}}(AZ)$, it just happens to be satisfied by the maximizer !
## 6.4 Weighted optimal projected explained variance

In order to select the SVD solution $Z = V_m$ among the maximizers of $expvar_{proj}^{opt}(AZ)$, we introduce weights $\mu_j$ such that:

$$ \mu_1 \ge \mu_2 \ge \dots \mu_m > 0, \quad (44) $$

and define a weighted optimal projected variance by:

$$ \expvar_{proj,\mu}^{opt}(AZ) = \max_{X^T X = I_m} \sum_{j=1 \dots m} \mu_j^2 \langle Az_j, x_j \rangle^2, $$

which coincides with $expvar_{proj}^{opt}(AZ)$ when $\mu_j = 1$ for all $j$. This leads to the weighted optimal projected explained variance block PCA formulation:

$$ \max_{\substack{\|z_j\| = 1 \\ j = 1 \dots m}} \expvar_{proj,\mu}^{opt}(AZ) = \max_{\substack{\|z_j\| = 1 \\ j = 1 \dots m}} \max_{X^T X = I_m} \sum_{j=1 \dots m} \mu_j^2 \langle Az_j, x_j \rangle^2 = \sum_{j=1 \dots m} \mu_j^2 \sigma_j^2. \quad (45) $$

The nice properties of this formulation are recalled in the next proposition:

**Proposition 6.3.** Let the singular values of A satisfy:

$$ \sigma_1 > \sigma_2 > \dots > \sigma_m > 0, $$

and the weights $\mu_j$ satisfy (44). Then the PCA loadings $Z = V_m$ and normalized components $X = U_m$ defined in (5) are one solution of the block PCA formulation (45) when the weights $\mu_j$ are constant, and its unique solution (up to a multiplication by $\pm 1$ of each column of course) when the weights $\mu_j$ are strictly decreasing, in which case the maximizers $Z^*$ and $X^*$ are independent of the weights $\mu_j$. The (unweighted) optimal projected variance explained by $Y^* = AZ^*$ is:

$$ \expvar_{proj}^{opt}(Y^*) = \sigma_1^2 + \dots + \sigma_m^2 \le \|A\|_F^2. $$

**Proof:** Exchanging the order of maximization in the center term of (45) solving analytically the maximization with respect to $Z$ gives:

$$ \max_{\substack{\|z_j\| = 1 \\ j = 1 \dots m}} \expvar_{proj,\mu}^{opt}(AZ) = \max_{X^T X = I_m} \sum_{j=1 \dots m} \mu_j^2 \|A^T x_j\|^2. \quad (46) $$

The last term in resp. (46) is the maximization of a weighted Rayleigh quotient for $A^T$, which is known to be equivalent to a PCA problem for $A^T$, and hence for $A$ (see for example Absil et al. [2008], recalled as Theorem 8.1 in the Appendix for the case of constant weights, and Brockett [1991] for the case of decreasing weights).

Of course, formulation (45) is of little interest for PCA itself, as there exists plenty of other efficient solution methods. But it will provide the starting point for the design of a sparse PCA algorithm, to be developed in the companion paper Chavent and Chavent [2023].
# 7 Conclusion

We have investigated the problem of defining the part of the variance of a data matrix explained by correlated components, such as those which arise when sparse loadings are searched for. We have established three compatibility conditions to be satisfied by any such explained variance definition in order to be compatible with the Principal Component Analysis (Condition 1), and to ensure a loss in explained variance when the components are correlated (Conditions 2 and 3). We have proved that the two existing and the four new definitions :

- Subspace (total variance of Shen and Huang [2008]),

- QR normalized,

- UP or Polar normalized,

- QR projected (adjusted variance of Zou et al. [2006]),

- UP or Polar projected,

- Optimal projected

all satisfy the two first compatibility conditions, but that only the three projected explained variance satisfy also the third one and provide proper explained variance definitions.

Numerical experiments have shown that the choice of a specific definition for the ranking of correlated components by explained variance is not critical. But we have shown that only the QR and the (weighted) optimal projected explained variance definitions admit the SVD solution as unique maximizer, and hence provide new explained variance block PCA formulations rid of orthogonality constraints on loadings. Their use for the construction of a group sparse PCA algorithm is the subject of a second paper Chavent and Chavent [2023].

# 8 Appendix

## 8.1 Generalized Rayleigh quotient

This is a classical result, see for example Absil et al. [2008] and Brockett [1991] :

**Theorem 8.1.** Let the loadings Z satisfy :

$$Z = [z_1 \dots z_m] \in \mathbb{R}^{p \times m}, \quad \mathrm{rank}(Z) = m \le \mathrm{rank}(A) \stackrel{\mathrm{def}}{=} r.$$

Then the generalized Rayleigh quotient

$$\mathrm{tr}\{(Z^T A^T A Z)(Z^T Z)^{-1}\}$$

satisfies :

$$\mathrm{tr}\{(Z^T A^T A Z)(Z^T Z)^{-1}\} \le \sigma_1^2 + \dots + \sigma_m^2 \le \|A\|_F^2,$$

and the left inequality becomes an equality if and only if :

$$\mathrm{span}\{Z\} = \mathrm{span}\{v_1 \dots v_m\},$$

where $v_1, \dots, v_m$ are the m first right singular vectors of A.
## 8.2 Proof of Lemma 3.1

Definition (12) of the subspace explained variance and the properties of the Rayleigh quotient $tr\{Z^T A^T AZ(Z^T Z)^{-1}\}$ recalled in Theorem 8.1 show that (9) and (10), and hence Properties 1 and 2, hold as well as (13) and (15).

It remains to prove (14) which shows that Condition 3 does not hold. So let $Y = AZ$ be orthogonal components :

$$ \langle y_j, y_k \rangle = 0, \quad j, k = 1 \dots m, j \neq k $$

corresponding to unit norm loadings :

$$ \|z_j\| = 1 \quad j = 1 \dots m, $$

and define $X, T$ by :

$$ x_j = y_j / \|y_j\|, \quad t_j = z_j / \|z_j\|, \quad j = 1 \dots m, $$

so that :

$$ X^T X = I_m. $$

Then on one side one has :

$$ \|Y\|_F^2 = \sum_{j=1 \dots m} \|y_j\|^2 = \sum_{j=1 \dots m} \frac{1}{\|t_j\|^2} = \mathrm{tr}\{\mathrm{diag}^{-1}(T^T T)\}, \quad (47) $$

and on the other side, as $Y$ and $X$ span the same subspace :

$$ \expvar_{subsp}(Y) = \expvar_{subsp}(X) = \mathrm{tr}\{(X^T X)(T^T T)^{-1}\} = \mathrm{tr}\{(T^T T)^{-1}\} \quad (48) $$

Formula (14) will be proved if we show that :

$$ \mathrm{tr}\{\mathrm{diag}^{-1}(T^T T)\} \le \mathrm{tr}\{(T^T T)^{-1}\}. \quad (49) $$

We use for that an idea taken from Miller [1969], and perform a QR-decomposition of $T$. By construction, the diagonal elements of $R$ satisfy :

$$ 0 < r_{i,i} \le \|t_i\|. $$

Then :

$$ \begin{aligned} T^T T &= R^T Q^T Q R = R^T R, \\ (T^T T)^{-1} &= R^{-1}(R^T)^{-1} = R^{-1}(R^{-1})^T, \end{aligned} $$

where $R^{-1}$ satisfies :

$$ R^{-1} = \text{upper triangular matrix, } [R^{-1}]_{i,i} = 1/r_{i,i}. $$

Hence the diagonal element of $(T^T T)^{-1}$ are given by : :

$$ \begin{aligned} [(T^T T)^{-1}]_{i,i} &= [R^{-1}(R^{-1})^T]_{i,i} \\ &= [R^{-1}]_{i,i}^2 + \sum_{j>i} [R^{-1}]_{i,j}^2 \\ &\ge [R^{-1}]_{i,i}^2 = 1/r_{i,i}^2 \ge 1/\|t_i\|^2. \end{aligned} \quad (50) $$
which gives (49) by summation over $i = 1...m$, and (14) is proved.

We suppose now that the orthogonal components $y_j, j = 1...m$ satisfy $\|Y\|_F^2 = \expvar_{subsp}(Y)$. Then (47) (48) imply that equality holds in (49) and hence all inequality in (50) are equalities :

1. first inequality : $[R^{-1}]_{i,j}^2 = 0$ for all $j > i \Rightarrow R^{-1}$ and hence $R$ are diagonal

2. second inequality : $1/r_{i,i}^2 = 1/\|t_i\|^2 \Rightarrow R$ is diagonal

But $R$ diagonal implies that the $t_j$ - and hence also the loadings $z_j$ - are orthogonal, which together with the hypothesis of orthogonal components $y_j$, implies that $(y_j/\|y_j\|, z_j)$ are pairs of singular vectors of $A$, which proves that $Z = V_m$ and ends the proof of (14).

## 8.3 Proof of Lemma 4.2

By construction, $\expvar_{proj}^X(AZ)$ satisfies clearly conditions 1 and 3 of Section 4. We prove now that it satisfies moreover (27), and hence also condition 2.

Let $\mathcal{E}^X$ be the ellipsoid of $\mathrm{span}(X) = \mathrm{span}(Y)$ image by $A$ of the unit sphere of $\mathrm{span}(Z)$. By construction one has :

$$y_j \in \mathcal{E}^X, \quad j = 1...m,$$

and the modified components $Y'$ defined by projection satisfy, c.f. (26) :

$$\|y'_{j}\| = |\langle y_j, x_j \rangle| \le \nu_j \stackrel{\mathrm{def}}{=} \max_{y \in \mathcal{E}^X} \langle y, x_j \rangle, \quad j = 1...m, \quad (51)$$

so that :

$$\expvar_{proj}^X Y \stackrel{\mathrm{def}}{=} \|Y'\|_F^2 \le \nu_1^2 + \dots + \nu_m^2. \quad (52)$$

We can now “box” the ellipsoid $\mathcal{E}^X$ in the parallelootope $\mathcal{P}^X$ of $\mathrm{span}\{X\}$ defined by :

$$\mathcal{P}^X = \{y \in \mathrm{span}\{X\} \mid -\nu_j \le \langle y, x_j \rangle \le +\nu_j, j = 1...m\},$$

(see figure 7). By construction, one can draw from each of the $2^m$ vertices of $\mathcal{P}^X$ $m$ orthogonal hyperplanes tangent to the ellipsoid $\mathcal{E}^X$, which implies that they are all on the orthoptic or Cartan sphere of the ellipsoid, whose radius is known to be the sum of the squares of the half principal axes $\sigma_j^X, j = 1...m$ of $\mathcal{E}^X$ (see for example the textbook Tauvel [2000]).

Hence :

$$\nu_1^2 + \dots + \nu_m^2 = (\sigma_1^X)^2 + \dots + (\sigma_m^X)^2. \quad (53)$$

Let then $y_1^X...y_m^X$ be vectors whose extremity are points of $\mathcal{E}^X$ located on its principal axes, so that :

$$\|y_j^X\| = \sigma_j^X, \quad j = 1...m, \quad \langle y_i^X, y_j^X \rangle = 0, \quad i, j = 1...m, i \ne j.$$

Property (14) of Lemma 3.1 applied to the orthogonal components $Y = Y^X$ gives :

$$(\sigma_1^X)^2 + \dots + (\sigma_m^X)^2 = \|Y^X\|^2 \le \expvar_{subsp} Y^X \le \sigma_1^2 + \dots + \sigma_m^2. \quad (54)$$

Combining inequalities (52) (53) (54) proves the inlem 4-2equality (27).
<figure><img src="image_9.png" /><figcaption>Figure 7: Illustration of the upper bound to $||Y'||_F^2$ in $\mathrm{span}\{Y\}$ when $Y'$ is defined by projection.</figcaption></figure>

## 8.4 Proof of Lemma 6.1 (notations of Section 8.3)

We prove first point 1 of the Lemma. Maximization of the convex function $X \rightsquigarrow \expvar_{proj} Y^X$ under the constraint $X^T X = I_m$ by algorithm (31) and passing to the limit proves (36).

We prove now the “only if” part of point 2. Let $X$ be given such that $\expvar_{proj} Y = \sigma_1^2 + \cdots + \sigma_m^2$. Then necessarily:

* equality holds in (54), and property (13) of subspace variance implies that the loadings $Z$ span the subspace $V_m$ of the $m$ first right singular vectors, which proves (37).

* equality holds in (51), which implies that for $j \neq k$ the normals to $\mathcal{E}^X$ at $y_j$ and $y_k$ are orthogonal (see Figure 7). The restriction $A_m$ of $A$ to $\mathrm{span}\{V_m\}$ is an isomorphism from $\mathrm{span}\{V_m\}$ to $\mathrm{span}\{U_m\}$, hence:

$$ \mathcal{E}^X = \{y \in \mathrm{span}\{U_m\} \mid \|A_m^{-1}y\|^2 = 1\}. $$

A normal $n(y)$ to $\mathcal{E}^X$ at $y$ is then:

$$ n(y) = \frac{1}{2} \nabla_y (\|A_m^{-1}y\|^2 - 1) = (A_m^{-1})^T A_m^{-1} y = (A_m^{-1})^T z, $$
and the orthogonality of $n(y_j)$ and $n(y_k)$ shows that :

$$ \langle n(y_j), n(y_k) \rangle = \langle (A_m^{-1})^T z_1, (A_m^{-1})^T z_2 \rangle = \langle z_1, (A_m^T A_m)^{-1} z_2 \rangle = 0, \quad (55) $$

which proves (38). When vectors and matrices are written on the singular bases $U_m$ and $V_m$, one has $(A_m^T A_m)^{-1} = \text{diag}\{\frac{1}{\sigma_1^2} \dots \frac{1}{\sigma_m^2}\}$.

We finally prove the “if” part of point 2. So let (37) (38) (39) hold. Property (37) implies that the half axes of $\mathcal{E}^X$ are $\sigma_1 \dots \sigma_m$, and (38) that the normal $n_j$ to $\mathcal{E}^X$ at $y_j, j = 1 \dots m$ are orthogonal. So one can box $\mathcal{E}^X$ with a parallelotope $\mathcal{P}^X$ with axes parallel to the normals $n_j$, and define $X$ as the orthonormal basis along its axes. Then the same reasoning as above for the proof of (27) shows that $\text{expvar}_{\text{proj}} Y = \sigma_1^2 + \dots + \sigma_m^2$, which ends the proof of the lemma.

## 8.5 Proof of Property (42)

Let the loadings $Z$ be such that the UP-projected explained variance of the components $Y = AZ$ satisfies $\text{expvar}_{\text{proj}}^{\text{UP}}(Y) = \sigma_1^2 + \dots + \sigma_m^2$. Then by definition of $\text{expvar}_{\text{proj}}^{\text{UP}}$ one has :

$$ \text{expvar}_{\text{proj}}^{\text{UP}}(Y) = \text{expvar}_{\text{proj}}^{X}(Y) \text{ with } Y = UP, U^T U = I_m, \quad P = P^T, P \ge 0 \quad (56) $$

But $\text{expvar}_{\text{proj}}^{X}(Y) = \sigma_1^2 + \dots + \sigma_m^2$, and point 1 of Lemma 6.1 implies the existence of $P'$ such that :

$$ Y \text{diag}(X^T Y) = XP' \text{ with } P'^T = P' \text{ and } P' \ge 0 \quad (57) $$

Comparison of the two last properties shows that :

$$ Y = UP = XP' \text{diag}(X^T Y)^{-1}, \quad (58) $$

and uniqueness of the polar decomposition implies that $P' \text{diag}(X^T Y)^{-1} = P$ and hence is symmetrical, which can happen in only two cases :

1. either $\text{diag}(X^T Y) = \lambda I_m$ for some $\lambda$, which gives the parasitic solution $Z = Z^\#$,

2. or $P'$ - and hence also $P$ itself - is diagonal, which implies that the components $Y = UP$ are orthogonal. But point 2 of Lemma 6.1 implies that they are also $(A_m^T A_m)^{-1}$ - orthogonal, which is possible only if the components $Y$ are proportional to the left singular basis $U_m$, which gives the SVD solution $Z \stackrel{\text{perm}}{=} V_m$.

## References

Pierre-Antoine Absil, Robert Mahony, and Rodolphe Sepulchre. *Optimization Algorithms on Matrix Manifolds*, volume 78. 12 2008. ISBN 978-0-691-13298-3. doi: 10.1515/9781400830244.
RW Brockett. Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems. *Linear algebra and its applications*, 146:79–91, 1991.

Guy Chavent. *Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications*. 01 2010. ISBN 978-90-481-2784-9. doi: 10.1007/978-90-481-2785-6.

Marie Chavent and Guy Chavent. A group sparse explained variance block pca. *Submitted*, 2023.

Michel Journée, Yurii Nesterov, Peter Richtárik, and Rodolphe Sepulchre. Generalized power method for sparse principal component analysis. *Journal of Machine Learning Research*, 11(Feb):517–553, 2010.

G Miller. Closed-form inversion of the gram matrix arising in certain least-squares problems. *IEEE Transactions on Circuit Theory*, 16(2):237–240, 1969.

R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL https://www.R-project.org/.

Haipeng Shen and Jianhua Z Huang. Sparse principal component analysis via regularized low rank matrix approximation. *Journal of multivariate analysis*, 99(6):1015–1034, 2008.

Patrice Tauvel. *Cours de géométrie: agrégation de mathématiques*. Dunod, Paris, 2000.

Hui Zou, Trevor Hastie, and Robert Tibshirani. Sparse principal component analysis. *Journal of computational and graphical statistics*, 15(2):265–286, 2006.