11institutetext: Michael Greenacre 22institutetext: Universitat Pompeu Fabra,
22email<EMAIL_ADDRESS>
# Compositional data analysis — linear algebra, visualization and
interpretation

Michael Greenacre

###### Abstract

Compositional data analysis is concerned with multivariate data that have a
constant sum, usually 1 or 100%. These are data often found in biochemistry
and geochemistry, but also in the social sciences, when relative values are of
interest rather than the raw values. Recent applications are in the area of
very high-dimensional “omics” data. Logratios are frequently used for this
type of data, i.e. the logarithms of ratios of the components of the data
vectors. These ratios raise interesting issues in matrix-vector
representation, computation and interpretation, which will be dealt with in
this chapter.

##  1 Introduction

Consider the table in Fig. 1(a): these are amounts spent in 2019 on four
different budget items in the European Union (EU), showing only six EU
countries with the addition of Iceland, Norway and Switzerland. The full data
set consists of 30 countries.

Figure 1: Expenditure in European Union countries plus Iceland, Norway and
Switzerland on four budget items, in 2019. In (a) the amounts are in millions
of local currency. In (b) the amounts are expressed as percentages relative to
the totals.

The amounts are in millions of local currency and are clearly not comparable
across the countries. Converting all the amounts to the same currency, for
example euros, alleviates but does not solve the data coding problem, since
some countries are small and others are large. Seeing that only relative
amounts are of principal interest, it seems that comparability across
countries is assured by simply expressing the four amounts in each row as
percentages, as in Fig. 1(b). The rows of data in Fig. 1(b) are called
compositions: nonnegative multivariate data with the constant sum constraint,
100% in this case. The act of transforming the original monetary values to
proportions, or percentages, is called in this context normalization or
closure.

But is the coding problem really solved? If additional budget items were added
to the table in Fig. 1(a), such as public services, social protection and
culture, the relative amounts in Fig. 1(b) would change, necessarily
decreasing, and they would reduce by different proportions since the
additional budget amounts would not be in proportion to the totals of the four
budget items shown in Fig. 1(a). This is what makes compositional data
different from any other multivariate data in Statistics — the values of each
component in the table depend on the values of the other components. It would
make no sense, for example, to compute correlations on such a data matrix,
since there are necessarily many negative correlations created by the constant
sum constraint, and the correlation between health and education would be
different in Fig. 1(b) from the correlation between health and education in an
expanded table of budget items also expressed as compositions.

Aitchison (1982, 1986) showed that using ratios of the components was a
solution to the data coding issue. Ratios remain constant between two
components irrespective of adding components to or removing components from a
composition — they are said to be subcompositionally coherent. Furthermore,
Aitchison proposed that ratios be logarithmically transformed: for example, if
$X_{j}$ and $X_{k}$ are two components then the logratio transformation is
$\log(X_{j}/X_{k})=\log(X_{j})-\log(X_{k})$, i.e. the difference in their
logarithms, which is a linear transformation on the log scale.

The objective of this chapter is to express the theory of logratio
transformations in linear algebra terms, and in the process give a flavour of
the visualization possibilities of logratios and their interpretation. In
Section 2, basic definitions and results are given in matrix-vector form, as
well as inverse logratio transformations, including the important topic of
log-contrasts. Section 3 shows how logratios can be visualized through cluster
analysis and biplots, and Section 4 concludes with a discussion.

##  2 Basic algebraic definitions and results

The practical aspects of compositional data analysis are given in two recent
books, by Greenacre (2018) and Filzmoser et al. (2018), and a comprehensive
review is given by Greenacre (2021). In this section, the algebra of
compositional data analysis is given in matrix–vector notation.

###  2.1 Logratio transformations and associated pattern matrices

Suppose ${\bf x}$ ($J\times 1$) is a $J$-component compositional vector of
positive values with sum 1: i.e. ${\bf
1}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu{\bf x}=1$, where
${\bf 1}$ is a vector of ones of appropriate order, in this case a vector of
$J$ ones. If $\log({\bf x})$ denotes the vector of log-transformed values,
then almost all logratio transformations are defined by a linear
transformation of the form ${\bf P}\log({\bf x})$, where ${\bf P}$ is called
the logratio pattern matrix. If ${\bf X}$ ($I\times J)$ denotes the data set,
with sampling units as rows and components as columns, then the compositions
are in the rows and the constant row sums are defined by postmultiplication by
${\bf 1}$: ${\bf X}{{\bf 1}}={{\bf 1}}$, where the ${\bf 1}$ on the right is
$J\times 1$. Similarly, if ${\bf L}=\log({\bf X})$ denotes the $I\times J$
matrix of logarithms of ${\bf X}$, then the application of the pattern matrix
to the rows of ${\bf L}$ implies post-multiplication on the right by the
transpose of the pattern matrix: ${\bf L}{\bf
P}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$. In the following
the matrix ${\bf P}$ will be subscripted by the type of logratio
transformation.

The pairwise logratio pattern matrix ${\bf P}_{\textsc{lr}}$
($\frac{1}{2}J(J-1)\times J$) is defined as:

${\bf P}_{\textsc{lr}}=\begin{bmatrix}[r]\ 1&-1&0&0&0&\cdots&0&\ \ 0&0\ \\\ \
1&0&-1&0&0&\cdots&0&\ \ 0&0\ \\\ \ 1&0&0&-1&0&\cdots&0&\ \ 0&0\ \\\ \
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\ \ \vdots&\vdots\ \\\ \
1&0&0&0&0&\cdots&0&\ \ 0&-1\ \\\ \ 0&1&-1&0&0&\cdots&0&\ \ 0&0\ \\\ \
0&1&0&-1&0&\cdots&0&\ \ 0&0\ \\\ \ 0&1&0&0&-1&\cdots&0&\ \ 0&0\ \\\ \
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\ \ \vdots&\vdots\ \\\ \
0&0&0&0&0&\cdots&1&\ \ 0&-1\ \\\ \ 0&0&0&0&0&\cdots&0&\ \ 1&-1\ \end{bmatrix}$
(1)

so that the $I\times\frac{1}{2}J(J-1)$ matrix of pairwise logratios is ${\bf
L}{\bf P}_{\textsc{lr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu$. Each row of ${\bf P}_{\textsc{lr}}$ defines a pairwise
logratio (LR), when applied to the logarithm of a compositio, $\log({\bf x})$.
Notice the lexicographic ordering in the rows, corresponding to the ratio
pairs (12), (13), (14), $\cdots$, (1$J$), (23), (24), (25), $\cdots$,
($J$–2,$\,J$), ($J$–1,$\,J$). We will abbreviate the term “pairwise logratio”
by LR throughout the rest of this chapter.

The matrix of LRs is of rank $J-1$, assuming $I\geq J$, otherwise it is of
rank $I-1$ (we shall assume for ease of description that there are at least as
many rows as columns in the compositional data matrix). This can be seen
easily using a result by Greenacre (2018, 2019) that a connected directed
acyclic graph (DAG) consisting of $J-1$ LRs generates all the
$\frac{1}{2}J(J-1)$ pairwise logratios through linear combinations. For
example, for the four-component example of Fig. 1(a) with six LRs, a possible
connected DAG is shown by the three solid arrows in Fig. 2, representing the
ratios Health/Economy, Health/Education and Defence/Education (the arrow
points towards the numerator component). The other three logratios, indicated
by dashed arrows, can be obtained from those in the DAG as ratios combined
either through multiplication or division depending on the direction of the
arrows: Education/Economy = (Health/Economy) / (Health/Education)), or in
linearized logratio form using addition or subtraction: log(Education) –
log(Economy = [log(Health)–log(Economy)] – [log(Health)–log(Education)].

-0.4cm

Figure 2: Directed acyclic graph (DAG) that connects the four components of
Fig. 1(a), indicated by the solid arrows. The dashed arrows indicate the other
three pairwise logratios that can be obtained from those of the DAG.

Assuming the four components above are in the order $\\{$Defence, Economy,
Health, Education$\\}$ as in the table in Fig. 1, the pattern matrix
associated with the three logratios (solid arrows) in the DAG above has this
form:

${\bf P}=\begin{bmatrix}[r]\ 0&-1&\ \ 1&0\ \\\ \ 0&0&\ \ 1&-1\ \\\ \ 1&0&\ \
0&-1\ \end{bmatrix}$ (2)

A special case of the LRs is the additive logratio (ALR) transformation, the
set of $J-1$ LRs where the denominator component is the same, called the
reference component. For example ,if the last component is the reference, then
the ALR pattern matrix ${\bf P}_{\textsc{alr}}$ ($(J-1)\times J$) is:

${\bf P}_{\textsc{alr}}=\begin{bmatrix}[r]\ 1&\ 0&\ 0&\cdots&\ 0&\ -1\ \\\ \
0&\ 1&\ 0&\cdots&\ 0&\ -1\ \\\ \ 0&\ 0&\ 1&\cdots&\ 0&\ -1\ \\\ \ \vdots&\
\vdots&\ \vdots&\ddots&\ \vdots&\ \vdots\ \\\ \ 0&\ 0&\ 0&\cdots&\ 1&\ -1\
\end{bmatrix}$ (3)

The matrix of ALRs, ${\bf L}{\bf
P}_{\textsc{alr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$
$(I\times(J-1))$, is of rank $J-1$. For example, if for the four components in
Fig. 1 the third component Health was chosen as the reference part, then (3)
would be a $3\times 4$ matrix with the $-1$s down the third column and a $1$
in each row in columns 1, 2 and 4. Fig. 3 shows the DAG associated with this
ALR transformation, where Health is placed in the centre and arrows emanate
outwards to the other three components.

-0.4cm

Figure 3: DAG corresponding to the ALR transformation of the four components
in Fig. 1, where Health is the reference part.

The next most important logratio transformation is the centered logratio (CLR)
transformation, the ratio of each component divided by the geometric mean of
all the components. The usual unweighted definition is the following, for a
row $[\,x_{1},x_{2},\ldots,x_{J}\,]$ of ${\bf X}$:

${\rm
CLR}(j)=\log\Biggl{(}\frac{x_{j}}{{\bigl{(}}\prod_{k}x_{k}{\bigr{)}}^{1/J}}\Biggr{)}=\log(x_{j})-\frac{1}{J}\sum_{k}\log(x_{k})\quad
j=1,\ldots,J$ (4)

but it is preferred here to give a more general weighted definition assuming
positive weights $c_{j}$ ($j=1,\ldots,J$) for the components, where
$\sum_{j}c_{j}=1$, and thus a weighted geometric mean in the denominator:

${\rm
CLR}(j)=\log\Bigl{(}\frac{x_{j}}{\prod_{k}x_{k}^{c_{k}}}\Bigr{)}=\log(x_{j})-\
\sum_{k}c_{k}\log(x_{k})\quad j=1,\ldots,J$ (5)

Hence, (4) is the special case of (5) with equal weights $1/J$ for all the $J$
components.

The CLR pattern matrix ${\bf P}_{\textsc{clr}}$ ($J\times J$) for the general
case is:

${\bf P}_{\textsc{clr}}=\begin{bmatrix}1-c_{1}&-c_{2}&-c_{3}&\cdots&-c_{J}\
\\\ -c_{1}&1-c_{2}&-c_{3}&\cdots&-c_{J}\ \\\
-c_{1}&-c_{2}&1-c_{3}&\cdots&-c_{J}\ \\\ \vdots&\vdots&\vdots&\ddots&\vdots\
\\\ -c_{1}&-c_{2}&-c_{3}&\cdots&1-c_{J}\ \end{bmatrix}$ (6)

Notice that ${\bf P}_{\textsc{clr}}$ is just the idempotent centering matrix
${\bf I}-{\bf 1}{\bf c}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu$, where ${\bf c}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu=[\,c_{1}\ c_{2}\cdots c_{J}\,]$. The rows of ${\bf
P}_{\textsc{clr}}$ sum to 0, and the matrix of CLRs, ${\bf L}{\bf
P}_{\textsc{clr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$,
has rank $J-1$, just like the matrix of LRs, ${\bf L}{\bf P}_{\textsc{lr}}$
and the matrix of ALRs, ${\bf L}{\bf P}_{\textsc{alr}}$.

More complex logratio transformations are the isometric logratios (ILRs)
(Egozcue and Pawlowsky-Glahn, 2005) and their slightly simpler special case,
the pivot logratios (PLRs) (Filzmoser et al., 2018) — these are often called
“balances” although the term can be misleading (Greenacre, 2020; Greenacre et
al., 2020). ILRs are log-transformed ratios of geometric means of groups of
components. In PLRs one of these groups is a single component and a linear
independent set of $J-1$ PLRs takes the components in a fixed order and the
numerator of each ratio is a single component and the denominator is the
geometric mean of the others “to the right of” the numerator component. For
both ILRs and PLRs there is a scalar constant involved, which is omitted here
for simplicity — see Greenacre (2018) for the exact (unweighted) definition.
Again it is preferred to give the more general weighted definition here, which
will be useful in practice when parts are considered differentially weighted.
The pattern matrix ${\bf P}_{\textsc{plr}}$ ($(J-1)\times J)$ for the weighted
case is:

${\bf P}_{\textsc{plr}}=\begin{bmatrix}\ 1&\ \ \
-\frac{c_{2}}{c_{2}+\cdots+c_{J}}&\ \
-\frac{c_{3}}{c_{2}+\cdots+c_{J}}&\cdots&-\frac{c_{J-2}}{c_{2}+\cdots+c_{J}}&\
\ -\frac{c_{J-1}}{c_{2}+\cdots+c_{J}}&\ \
-\frac{c_{J}}{c_{2}+\cdots+c_{J}}\\\\[5.16663pt] \ 0&\ \ \ 1&\ \
-\frac{c_{3}}{c_{3}+\cdots+c_{J}}&\cdots&-\frac{c_{J-2}}{c_{3}+\cdots+c_{J}}&\
\ -\frac{c_{J-1}}{c_{3}+\cdots+c_{J}}&\ \
-\frac{c_{J}}{c_{3}+\cdots+c_{J}}\\\\[5.16663pt] \ 0&\ \ \ 0&\ \
1&\cdots&-\frac{c_{J-2}}{c_{4}+\cdots+c_{J}}&\ \
-\frac{c_{J-1}}{c_{4}+\cdots+c_{J}}&\ \
-\frac{c_{J}}{c_{4}+\cdots+c_{J}}\\\\[5.16663pt] \vdots&\ \ \
\vdots&\vdots&\ddots&\vdots&\ \ \vdots&\ \ \vdots\\\\[5.16663pt] \ 0&\ \ \ 0&\
\ 1&\cdots&1&-\frac{c_{J-1}}{c_{J-1}+c_{J}}&\ \
-\frac{c_{J}}{c_{J-1}+c_{J}}\\\\[5.16663pt] \ 0&\ \ \ 0&\ \ 0&\cdots&0&1&\ \
-1\end{bmatrix}$ (7)

Like the pattern matrices before, (1), (3) and (6), the row sums of ${\bf
P}_{\textsc{plr}}$ are zero and the matrix ${\bf L}{\bf
P}_{\textsc{plr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ of
PLRs has rank $J-1$. There are $J!$ permutations of the $J$ components, hence
$J!$ sets of PLRs possible, depending on the ordering of the $J$ components.
For the four components of Fig. 1, in the order given, a graph of the
associated PLRs is in the form of the binary dendrogram of Fig. 4. Since there
are $4!=24$ ordered permutations possible for this example, there are 24
different sets of PLRs possible. Notice that the last member of a set of PLRs
is a simple LR.

-0.4cm

Figure 4: Dendrogram graph associated with three PLRs corresponding to the
order of the components Defence, Economy, Health, Education.

The ILRs involving more general ratios of geometric means is usually defined
according to a dendrogram graph, and the number of possible dendrograms
increases even more rapidly with the number of components. The ILR pattern
matrix becomes more difficult to express in general, so we will just give the
special case associated with the dendrogram in Fig. 5, assuming the order
shown of the components.

-0.4cm

Figure 5: Dendrogram graph associated with three ILRs corresponding to the
three ratios (moving down the tree): (1) geometric mean of Defence and Economy
/ geometric mean of Health and Education, (2) Defence / Economy, (3) Health /
Education.

The specific pattern $3\times 4$ matrix for this four-component example
corresponding to Fig. 5, again with rows summing to 1, would thus be:

${\bf P}_{\sf ILR}=\begin{bmatrix}\frac{c_{1}}{c_{1}+c_{2}}&\ \
\frac{c_{2}}{c_{1}+c_{2}}&\ \ -\frac{c_{3}}{c_{3}+c_{4}}&\ \
-\frac{c_{4}}{c_{3}+c_{4}}\\\\[3.01385pt] 1&\ \ -1&\ \ \ \ 0&\ \ \ \
0\\\\[3.01385pt] 0&\ \ \ \ 0&\ \ \ \ 1&\ \ -1\end{bmatrix}$ (8)

Finally, there is a class of nonlinear transformations called amalgamation (or
summated) logratios, abbreviated as SLRs. Like the LRs, these are true
balances between the components as they use sums rather than geometric means
when combining components in the ratios (Greenacre, 2020). For example, the
dendrogram in Fig. 5 would translate to the following three SLRs, the last two
of which are regular LRs:

$(1)\log\left(\frac{\textrm{Defence+Economy}}{\textrm{Health+Education}}\right)\quad(2)\log\left(\frac{\textrm{Defence}}{\textrm{Economy}}\right)\quad(3)\log\left(\frac{\textrm{Health}}{\textrm{Education}}\right)$
(9)

The first SLR above cannot be written as a linear function of the logarithms,
but has the advantage of being more easily interpreted in a practical
application. It is still isomorphic, however, as shown in the next subsection.

###  2.2 Inverting logratio transformations

Each of the logratio transformations defined in the previous subsection can be
inverted back to the original compositions, including the nonlinear SLRs. Each
in turn relies on a square inversion pattern matrix, denoted by ${\bf Q}$,
closely related to the respective logratio pattern matrix ${\bf P}$. The back-
transformation in each case involves solving a system of $J$ linear equations
of the form ${\bf Q}{\bf x}={\bf e}$, i.e. ${\bf x}={\bf Q}^{-1}{\bf e}$,
where ${\bf Q}$ involves elements that are functions of the logratios, and
${\bf e}$ is the vector $[\,0\ \,0\,\cdots\,0\
\,1\,]^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$.

Starting with the ALRs, suppose that ${\bf y}={\bf P}_{\textsc{alr}}{\bf x}$
is the vector of $J-1$ ALRs, where ${\bf x}$ is a $J\times 1$ composition and
${\bf P}_{\textsc{alr}}$ is given by (3). Then the inverse operation of
finding ${\bf x}$ from ${\bf y}$ is the solution of the following equation:

${{\bf Q}}_{\textsc{alr}}{\bf x}=\begin{bmatrix}\ 1&\ 0&\cdots&0&\
-e^{y_{1}}\\\ \vspace{-0.1cm}\ 0&\ 1&\cdots&0&\ -e^{y_{2}}\\\ \ \vdots&\
\vdots&\ddots&\vdots&\vdots\\\ \ 0&\ 0&\cdots&1&\ -e^{y_{J-1}}\\\ \ 1&\
1&\cdots&1&\ 1\end{bmatrix}{\bf x}=\begin{bmatrix}0\\\ \vspace{-0.1cm}0\\\
\vdots\\\ 0\\\ 1\end{bmatrix}$ (10)

Notice that the $-1$s in the pattern matrix, i.e. in the denominator
positions, are substituted by $-e^{y_{1}},-e^{y_{2}},\ldots,-e^{y_{J-1}}$ in
rows $1,2,\ldots,J-1$ respectively, and then a row of $1$s is added, which
with the last $1$ in the right hand vector imposes the constraint that the
$x_{j}$s in the solution sum to 1 (Greenacre, 2020). This is equivalent to the
simpler calculation of exponentiating the $J-1$ logratios and expanding them
with a 1, i.e. $[\ e^{y_{1}}\ e^{y_{2}}\cdots\ e^{y_{J-1}}\ \ 1\,]$ and then
normalizing the result to sum to 1. So, although using the linear equations
approach by solving (10) is actually an inefficient way of back-transforming
the compositions, it is enlightening because it provides a way to deal with
any set of $J-1$ independent LRs, e.g. the three LRs defined in the DAG of
Fig. 2.

The way of inverting a set of $J-1$ independent pairwise logratios (LRs) is a
simple generalization of the system of equations (10) above. The $(J-1)\times
J$ pattern matrix ${\bf P}_{\textsc{lr}}$ for the transformation to the LRs,
now has each of the $J-1$ rows corresponding to a specific LR,
$\log(x_{j}/x_{j^{\prime}})$, with a $1$ in column $j$ of the numerator part
and $-1$ in column $j^{\prime}$ of the denominator part. To obtain the
inversion pattern matrix, the $-1$ is again replaced in row $k$ by
$-e^{y_{k}}$ and then a row of $1s$ is added as the last row, as before. For
example, for the DAG in Fig. 2, the logratio pattern matrix and corresponding
inversion pattern matrix, with their matrix equations are, in the order
Defence, Economy, Health, Education, and LRs Health/Economy, Health/Education
and Defence/Education, for ${\bf x}$ ($4\times 1$) a compositional vector,
${\bf y}$ ($3\times 1$) the logratio transformation, and ${\bf e}=[\,0\,\ 0\,\
0\,\ 1\,]^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$:

${\bf P}_{\textsc{lr}}=\begin{bmatrix}[r]0&-1&\ 1&0\\\ 0&0&\ 1&-1\\\ 1&0&\
0&-1\end{bmatrix},\ \ {\rm i.e.\ }{\bf y}={\bf P}_{\textsc{lr}}\log({\bf
x})\quad{\bf Q}_{\textsc{lr}}=\begin{bmatrix}0&-e^{y_{1}}&1&0\\\
0&0&1&-e^{y_{2}}\\\ 1&0&0&-e^{y_{3}}\\\ 1&1&1&1\end{bmatrix},\ \ {\rm i.e.\
}{\bf x}={\bf Q}_{\textsc{lr}}^{-1}{\bf e}$ (11)

The same approach can be used to invert a set of $J-1$ independent
amalgamation logratio balances (SLRs), which must involve each part at least
once. The SLRs have a pattern matrix ${\bf P}_{\textsc{slr}}$ indicating which
parts are in the numerator and in the denominator, but this is not the matrix
of the transformation, which is not linear in $\log({\bf x})$, hence not a
contrast matrix. For example, the pattern matrix for the SLRs in (9) as well
as the corresponding inversion pattern matrix are as follows:

${\bf P}_{\textsc{slr}}=\begin{bmatrix}[r]1&1&-1&-1\\\ 1&-1&0&0\\\
0&0&1&-1\end{bmatrix}\qquad{\bf
Q}_{\textsc{slr}}=\begin{bmatrix}1&1&-e^{y_{1}}&-e^{y_{1}}\\\
1&-e^{y_{2}}&0&0\\\ 0&0&1&-e^{y_{3}}\\\ 1&1&1&1\end{bmatrix}$ (12)

Each row of ${\bf P}_{\textsc{slr}}$ corresponding to an SLR has $1$s in the
columns of the numerator parts, and $-1$s in the columns of the denominator
parts. Then, as before, replace all the $-1$s in the $k$-th row of ${\bf
Q}_{\textsc{slr}}$ by $-e^{y_{k}}$ and add a row of $1s$ as the last row,
which again automatically closes the parts in the solution. The compositional
vector can then be recovered by ${\bf x}={\bf Q}_{\textsc{lsr}}^{-1}{\bf e}$.
Notice that both the LR and ALR inverse transformations are special cases of
the SLR one, where LRs and ALRs have only one numerator and one denominator
part. For specific examples of inverse transforms of sets of LRs and sets of
SLRs, see Greenacre (2020).

###  2.3 Log-contrasts

A log-contrast is a linear combination of logarithms of all the components of
composition, with the condition that the coefficients sum to $0$:

$\sum_{j}a_{j}\log(x_{j}),\quad{\rm where\ }\sum_{j}a_{j}=0.$ (13)

The logratio pattern matrices ${\bf P}_{\textsc{alr}}$ in (3), ${\bf
P}_{\textsc{clr}}$ in (6), and ${\bf P}_{\textsc{plr}}$ in (7), as well as a
particular pattern matrix associated with a DAG such as (2), transform the
vector $\log({\bf x})$ to the corresponding set of $J-1$ logratios, or $J$
logratios in the case of the CLR transformation. For any one of these sets of
logratios, denoted in general by $\ell_{1},\ell_{2},\ldots$, the coefficients
$\bf c$ of a linear combination of them $c_{1}\ell_{1}+c_{2}\ell_{2}+\cdots$
can be converted to the coefficients of the log-contrast simply by pre-
multiplying $\bf c$ by the transpose of the pattern matrix. For example, for a
linear combination of ALRs, the $J$ coefficients of the log-contrast are ${\bf
a}={\bf P}_{\textsc{alr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu{\bf c}$. This result is useful when a linear combination of
logratios, used as explanatory variables in a generalized linear model, is
estimated in explaining/predicting a response variable (Coenders and
Pawlowsky-Glahn, 2020).

This result is illustrated for four different transformations, for the linear
modelling of a response variable in the form of the proportion of total budget
in the same 30 countries devoted to Housing and Community Amenities,
logarithmically transformed, denoted by $\log(y)$. For example, if one defines
the ALR with Health as the reference component, then the estimated regression
model is, along with p-values in parentheses and proportion of explained
variance, $R^{2}$:

$\log(y)=-0.482+0.029\log\left(\frac{\textrm{\scriptsize
Defence}}{\textrm{\scriptsize Health
}}\right)+1.032\log\left(\frac{\textrm{\scriptsize
Economy}}{\textrm{\scriptsize
Health}}\right)-0.904\log\left(\frac{\textrm{\scriptsize
Education}}{\textrm{\scriptsize Health}}\right)$ $\hskip
54.06006pt(p=0.004)\hskip 34.14322pt(p=0.030)\hskip 34.14322pt(p=0.83)\hskip
28.45274ptR^{2}=0.295$

The coefficients of the log-contrast, using the corresponding ALR pattern
matrix, are then (always remembering the order of the components):

${\bf a}={\bf P}_{\textsc{alr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu{\bf c}=\begin{bmatrix}[r]1&0&0\ \\\ 0&1&0\ \\\ -1&-1&-1\ \\\
0&0&1\ \end{bmatrix}\begin{bmatrix}[r]0.029\ \\\ 1.032\ \\\ -0.904\
\end{bmatrix}=\begin{bmatrix}[r]0.029\ \\\ 1.032\ \\\ -0.157\ \\\ -0.904\
\end{bmatrix}$

It can be verified that ${\bf 1}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu{\bf a}=1$. Thus the regression model can be written as the
constant plus the log-contrast:

$\log(y)=-0.482+0.029\log(\textrm{Defence})+1.032\log(\textrm{Economy})-0.157\log(\textrm{Health})-0.904\log(\textrm{Education})$

From the form of ${\bf
P}_{\textsc{alr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ and
the fact that any ALR transformation will give the same log-contrast, it is
clear that the coefficient of the reference part is the one that will change.
For example, if Education is the reference, the regression coefficients turn
out to be $[\ 0.029\ \ 1.032\ \ -0.157\
]^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ and the log-
contrast is shown to be identical:

${\bf a}={\bf P}_{\textsc{alr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu{\bf c}=\begin{bmatrix}[r]1&0&0\ \\\ 0&1&0\ \\\ 0&0&1\ \\\
-1&-1&-1\ \end{bmatrix}\begin{bmatrix}[r]0.029\ \\\ 1.032\ \\\ -0.157\
\end{bmatrix}=\begin{bmatrix}[r]0.029\ \\\ 1.032\ \\\ -0.157\ \\\ -0.904\
\end{bmatrix}$

For any set of pairwise logratios, once again the same log-contrast is
obtained. For example, here is the regression model using the logratios in the
DAG of Fig. 2:

$\log(y)=-0.482-1.032\log\left(\frac{\textrm{\scriptsize
Health}}{\textrm{\scriptsize Economy
}}\right)+0.875\log\left(\frac{\textrm{\scriptsize
Health}}{\textrm{\scriptsize
Education}}\right)+0.029\log\left(\frac{\textrm{\scriptsize
Defence}}{\textrm{\scriptsize Education}}\right)$ $\hskip
59.75095pt(p=0.004)\hskip 34.14322pt(p=0.036)\hskip 39.83368pt(p=0.83)\hskip
22.76228ptR^{2}=0.295$

Then, using the pattern matrix ${\bf P}$ in (2):

${\bf a}={\bf P}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu{\bf
c}=\begin{bmatrix}[r]0&0&1\ \\\ -1&0&0\ \\\ 1&1&0\ \\\ 0&-1&-1\
\end{bmatrix}\begin{bmatrix}[r]-1.032\ \\\ 0.875\ \\\ 0.029\
\end{bmatrix}=\begin{bmatrix}[r]0.029\ \\\ 1.032\ \\\ -0.157\ \\\ -0.904\
\end{bmatrix}$

The same result is obtained for the CLRs, as well as any set of ILRs or PLRs.
One difference with the CLRs is that the pattern matrix is $4\times 4$ and
only three coefficients are obtained in a regression, so the fourth one has to
be set to zero. In all cases the constant as well as the $R^{2}$ and the
p-value for the whole model (which is $p=0.026$) are identical across the
variations, as they all reduce to the same log-contrast.

Often, some type of variable selection is made to arrive at a more
parsimonious model. Selecting fewer explanatory variables implies forming a
subcomposition of the parts. A statistical criterion is needed to make the
selection and there are many possible ways to do achieve this. For example,
one could do a permutation test on the coefficients of the log-contrast. Using
the CLR transformation, and randomizing the order of the response variable 999
times, the p-values for each log-contrast coefficient were estimated as:

Defence: $p=0.86$, Economy: $p=0.003$, Health: $p=0.66$, Education: $p=0.037$

.

It seems that only one logratio, that of Economy/Education, can be used as a
predictor, and gives the following result:

$\log(y)=-0.562+1.010\log\left(\frac{\textrm{\scriptsize
Economy}}{\textrm{\scriptsize Eduction }}\right)$ $\hskip
56.9055pt(p=0.002)\hskip 34.14322ptR^{2}=0.286$

This leads to the trivial log-contrast of the two components in the model:

$\log(y)=-0.562+1.010\log(\textrm{Economy})-1.010\log(\textrm{Education})$

For more details about logratios used as predictors in linear modelling, see
Coenders and Pawlowsky-Glahn (2020).

##  3 Logratio visualization

In this section we look at various ways of visualizing a compositional data
set. Basically, once a logratio transformation is made, any of the various
well-known multivariate visualization methods can be implemented, such as
cluster analysis and dimension-reduced component methods. Care has to be taken
in the interpretation because of the unit sum constraint on the original data.
Since these methods rely on interpoint distances, the first thing to do is to
define the distance measures between rows and between columns of the data
matrix.

If the matrix ${\bf Z}=[z_{i,jj^{\prime}}]$ ($I\times\frac{1}{2}J(J-1)$)
denotes the matrix of LRs $\log(x_{ij}/x_{ij^{\prime}})$, and ${\bf
Y}=[y_{ij}]$ ($I\times J$) the CLR-transformed data set, then the _logratio
distance_ $d_{ii^{\prime}}$ between samples $i$ and $i^{\prime}$, can be
defined in two equivalent forms, shown for the weighted and unweighted
versions in (14) and (15) respectively (Greenacre, 2018):

$d_{ii^{\prime}}=\sqrt{{\sum\sum}_{j<j^{\prime}}c_{j}c_{j^{\prime}}(z_{i,jj^{\prime}}-z_{i^{\prime},jj^{\prime}})^{2}}=\sqrt{\sum_{j}c_{j}(y_{ij}-y_{ij^{\prime}})^{2}}$
(14)
$d_{ii^{\prime}}=\sqrt{\frac{1}{J^{2}}{\sum\sum}_{j<j^{\prime}}(z_{i,jj^{\prime}}-z_{i^{\prime},jj^{\prime}})^{2}}=\sqrt{\frac{1}{J}{\sum}_{j}(y_{ij}-y_{ij^{\prime}})^{2}}$
(15)

The advantage of the versions using CLRs is the use of a much narrower matrix,
but — as will be emphasized repeatedly — the results should always be
interpreted in terms of pairwise logratios. The CLRs as such have no inherent
interpretation as representing the components and simply act as a short cut to
analysing all the pairwise logratios.

There are similar results for the columns, by transposing the data set,
renormalizing and performing the same operations. The samples are usually
unweighted (i.e., with weights $1/I$) but can also be differentially weighted
if required — see Greenacre (2018). Later for the definition of the biplot the
theory is presented in complete generality with weights on the rows and the
columns.

Fig. 7 shows the clustering the rows and columns respectively of an extended
data set of country budget items (Fig. 6), using the easyCODA package
(Greenacre, 2018) in R (R Core Team, 2021), and Ward clustering in each case
(Ward, 1963).

Figure 6: Expenditure in European Union countries plus Iceland, Norway and
Switzerland on ten budget items, in 2019. The data are expressed as
percentages of the expenditure on these items (i.e., row sums are 100%). Some
columns names have been slightly abbreviated — see the original longer names
in Fig. 7 below. Figure 7: Clustering of (a) countries and (b) budget items,
using Ward clustering of the logratio distances. The budget items are
weighted, as in (14), using weights proportional to their average marginal
proportions.

A completely different way of performing the clustering is by successively
amalgamating the components (Greenacre, 2020), shown in Fig. 8. Both
unweighted and weighted versions are shown, giving different results, and
different from the Ward clustering in Fig. 7(b). The unweighted version is
identical to the graph structure of a set of PLRs, while the weighted version
de-emphasizes the role of Defence, and shows that there are three pairs of
items that could be merged: Housing and community amenities with Economic
affairs, Environmental protection with Public order and safety, and Education
with Recreation, culture and religion.

Figure 8: Clustering of budget items using amalgamation clustering: (a)
unweighted and (b) weighted.

Next, a logratio biplot is shown in Fig. 9, using logratio analysis (LRA).
This biplot is obtained using the singular value decomposition (SVD) of the
double-centered matrix of the log-transformed compositional data set,
$\log({\bf X})$. An unweighted or weighted version of LRA is possible, and
here weighted LRA is used, where weights are imposed on the budget items equal
to the average proportions of the items across the countries (Greenacre and
Lewi, 2009; Greenacre, 2018). This weighting will reduce the influence of some
budget items with low average proportions but high logratio variance. The
configuration of the countries is approximating the weighted logratio
distances in (14).

Figure 9: Weighted logratio biplot of the data in Fig. 6, with asymmetric
biplot scaling: rows in principal coordinates, columns in standard
coordinates.

Using matrix-vector notation, the sequence of steps to perform LRA and arrive
at the biplot in Fig. 9 is as follows, starting from the matrix of log-
transformed compositional data, $\log({\bf X})$. Note that the most general
row- and column-weighted version is given here, with row and column weights
${\bf r}$ and ${\bf c}$ respectively. Usually, but not necessarily, the rows
are equally weighted: ${\bf r}=(1/I){\bf 1}$, and the description weighted or
unweighted LRA refers to the nature of the column weights ${\bf c}$, which can
be different (e.g., by default equal to the marginal average proportions of
the components (Greenacre and Lewi, 2009)) or equal: ${\bf c}=(1/J){\bf 1}$
(Aitchison and Greenacre, 2002).

$\displaystyle\textrm{Double-centre the matrix}\log({\bf X}):$
$\displaystyle\mkern-1.0mu{\bf Z}\mkern-2.0mu=({\bf I}-{\bf
1r}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu)\log({\bf X})({\bf
I}-{\bf 1c}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu)^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ (16)
Apply weights to rows and columns: $\displaystyle{\bf S}={\bf
D}_{r}^{\frac{1}{2}}{\bf Z}{\bf D}_{c}^{\frac{1}{2}}$ (17) Perform the SVD:
$\displaystyle{\bf S}={\bf UD}_{\alpha}{\bf
V}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ (18) Principal
coordinates of rows: $\displaystyle{\bf F}={\bf D}_{r}^{\\!-\frac{1}{2}}{\bf
UD}_{\alpha}$ (19) Standard coordinates of columns:
$\displaystyle{\mathbf{\Gamma}}={\bf D}_{c}^{\\!-\frac{1}{2}}{\bf V}$ (20)
Contribution coordinates of columns:
$\displaystyle{\mathbf{\Gamma}}^{\mkern-1.0mu\ast}\mkern-7.0mu={\bf
D}_{c}^{\frac{1}{2}}{\mathbf{\Gamma}}={\bf V}$ (21)

Two-dimensional biplots of the results are given by plotting the row principal
coordinates in the first two columns of $\bf F$ with either the corresponding
column standard coordinates in the first two columns of ${\mathbf{\Gamma}}$
(asymmetric biplot) or those of the contribution coordinates in
${\mathbf{\Gamma}}^{\mkern-1.0mu\ast}$ (contribution biplot).

The two-dimensional symmetric map of the LRA is the plotting of the first two
columns of the row principal coordinates F in (19) jointly with those of the
column principal coordinates G:

Principal coordinates of columns: $\displaystyle{\bf G}={\bf
D}_{c}^{\\!-\frac{1}{2}}{\bf VD}_{\alpha}$ (22)

This is not a biplot, strictly speaking, but has the practical advantage that
the row and column points are equally scaled along the principal axes, their
weighted variances both being equal to the amount of variance explained on the
axes, i.e. the eigenvalues or squared singular values $\alpha$${}_{k}^{2}$ on
axis $k$, $k=1,2$.

The sum of the squared singular values, $\sum_{k}$$\alpha$${}_{k}^{2}$ (i.e.
sum of eigenvalues), equals the total logratio variance and
$\alpha$${}_{k}^{2}$ is the part of variance explained by axis $k$, usually
expressed as a percentage of this total, as shown on the axes in Fig. 9.

The steps (16)–(21) are equivalent to performing a principal component
analysis (PCA) on the matrix of (weighted) CLRs, because the CLRs are the row-
centered $\log({\bf X})$ (i.e., $\log({\bf X})({\bf I}-{\bf
1c}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu)^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ in
(16)) and PCA automatically performs column-centering (i.e., the centering
$({\bf I}-{\bf 1r}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu)$
in (16), with equal weights in ${\bf r}$), hence the double-centering. The
steps from (17) onwards define the PCA with its variations of display
coordinates. Because of the double-centering of $\log({\bf X})$, the
definition of each principal component as a linear combination of the CLRs
turns out to be a log-contrast.

The interpretation of Fig. 9 is not the same as a regular PCA, however. The
CLRs as variables are not interpretable per se, but rather the differences
between pairs of CLR points, which depict the pairwise logratios themselves.
Aitchison and Greenacre (2002) show that the LRA biplot optimizes the display
of these pairwise logratios, which is not the case in a regular PCA (i.e., in
a regular PCA the optimization of the variables is not equivalent to the
optimization of the differences between pairs of variables). Thus, to the
horizontal dispersion in Fig. 9 is due to logratios such as Public Order and
Safety divided by Social Protection, while the vertical dispersion is
dominated by the logratio of Defence relative to Health, for example, where it
is clear that Iceland’s ratios of Defence relative to the other budget items
are low.

The equivalence between the PCA of the CLRs and the PCA of the LRs (i.e., LRA
in both cases), can be neatly shown using the respective logratio pattern
matrices defined earlier. The proof follows the one given for unweighted
logratios in the Appendix of Aitchison and Greenacre (2002), but is more
elegantly defined for the general weighted case. The double-centred matrix
${\bf Z}$ in (16) is the matrix ${\bf L}=\log({\bf X})$ post-multiplied by the
transposed column-centering matrix, which is identical to the transposed CLR
pattern matrix in (6), and pre-multiplied by the row-centering matrix, i.e.
${\bf Z}=({\bf I}-{\bf 1}{\bf r}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu){\bf L}{\bf
P}_{\textsc{clr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$.
The corresponding result for the matrix of pairwise LRs ${\bf L}{\bf
P}_{\textsc{lr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$, is
the row-centered matrix ${\bf Y}=({\bf I}-{\bf 1}{\bf
r}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu){\bf L}{\bf
P}_{\textsc{lr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$. The
weights in the two respectve cases are $c_{1},c_{2},\ldots,c_{J}$ for the $J$
CLRs and $c_{1}c_{2},c_{1}c_{3},\ldots,c_{J-1}c_{J}$ for the
$\frac{1}{2}J(J-1)$ LRs, gathered in the diagonals of the diagonal matrices
${\bf D}_{c}$ and ${\bf D}_{cc}$ respectively. The matrix of weighted scalar
products for the rows of ${\bf Z}$ and ${\bf Y}$ are then identical. These
$I\times I$ matrices are called the form matrices:

${\bf Z}{\bf D}_{c}{\bf Z}^{\sf\scriptsize\raisebox{0.90417pt}{\sf
T}}\mkern-1.0mu={\bf Y}{\bf D}_{cc}{\bf
Y}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu$ (23)

thanks to the result ${\bf
P}_{\textsc{clr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu{\bf
D}_{c}{\bf P}_{\textsc{clr}}={\bf
P}_{\textsc{lr}}^{\sf\scriptsize\raisebox{0.90417pt}{\sf T}}\mkern-1.0mu{\bf
D}_{cc}{\bf P}_{\textsc{lr}}$. This proves that the configuration of the
samples is identical. The distances in (14) correspond exactly to the form
matrices in (23).

Notice that the equivalence of the PCA of the CLRs and the PCA of the LRs is a
particular result thanks to the definition of the CLRs. In regular PCA the
low-dimensional result is not optimal for differences between variables. For
example, Gabriel (1972) gives examples of difference vectors in biplots but
makes no statements about their optimality, whereas Greenacre (2003) shows
specifically how difference vectors can be optimally displayed.

##  4 Discussion

In compositional data analysis it is all the pairwise logratios that are
analysed and interpreted simultaneously. For multivariate analysis the
centered logratio transformation is the most useful and most efficient, since
it implies the analysis of all the pairwise logratios. The matrix of CLRs has
$J$ columns but is of rank $J-1$ (it is assumed that the number of rows $I$ is
greater than the number of columns $J$, otherwise the rank is $I-1$). When a
multivariate analysis such as discriminant analysis requires the inversion of
a singular covariance matrix, then the generalized inverse needs to be used.
When there is a response variable and the compositional variables are used as
explanatory variables (see Section 2.6), there are important issues of effect-
size interpretation, discussed by Coenders and Pawlowsky-Glahn (2020).

As John Aitchison said, “Compositional data analysis is simple” (Aitchison,
1997). The basic concept and initial step is to use the logratio
transformation, after which statistical analysis proceeds very much as before.
But care needs to be taken in the interpretation: the analyst has to think in
terms of pairwise logratios and realize the implications on the results of the
initial normalization of the data to have constant sums.

## References

  * Aitchison [1982] J. Aitchison.  The statistical analysis of compositional data (with discussion).  _J R Stat Soc Ser B_ , 44:139–77, 1982.
  * Aitchison [1986] J. Aitchison.  _The Statistical Analysis of Compositional Data_.  Chapman & Hall, London, 1986.
  * Aitchison [1997] J. Aitchison.  The one-hour course in compositional data analysis, or compositional data analysis is simple.  In V. Pawlowsky-Glahn, editor, _Proceedings of IAMG’97_ , pages 3–35. International Association for Mathematical Geology, 1997.
  * Aitchison and Greenacre [2002] J. Aitchison and M. Greenacre.  Biplots of compositional data.  _J R Stat Soc Ser C (Appl Stat)_ , 51:375–92, 2002.
  * Coenders and Pawlowsky-Glahn [2020] G. Coenders and V. Pawlowsky-Glahn.  On interpretations of tests and effect sizes in regression models with a compositional predictor.  _SORT_ , 44:201–20, 2020.
  * Egozcue and Pawlowsky-Glahn [2005] J.J. Egozcue and V. Pawlowsky-Glahn.  Groups of parts and their balances in compositional data analysis.  _Math Geol_ , 37:795–828, 2005.
  * Filzmoser et al. [2018] P. Filzmoser, K. Hron, and M. Templ.  _Applied Compositional Data Analysis_.  Oxford University Press, Oxford, 2018.
  * Gabriel [1972] K.R. Gabriel.  Analysis of meteorological data by means of canonical decomposition and biplots.  _J Appl Meteor Climat_ , 11:1077–7, 1972.
  * Greenacre [2003] M. Greenacre.  Singular value decomposition of matched matrices.  _J Appl Stat_ , 30(10):1101–13, 2003.
  * Greenacre [2018] M. Greenacre.  _Compositional Data Analysis in Practice_.  Chapman & Hall / CRC Press, Boca Raton, Florida, 2018.
  * Greenacre [2019] M. Greenacre.  Variable selection in compositional data analysis using pairwise logratios.  _Math Geosc_ , 51:649–82, 2019.
  * Greenacre [2020] M. Greenacre.  Amalgamations are valid in compositional data analysis, can be used in agglomerative clustering, and their logratios have an inverse transformation.  _Appl Comput Geosc_ , 5:100017, 2020.  URL https://doi.org/10.1016/j.acags.2019.100017.
  * Greenacre [2021] M. Greenacre.  Compositional data analysis.  _Annu Rev Stat Appl_ , 8:271–99, 2021.
  * Greenacre and Lewi [2009] M. Greenacre and P. Lewi.  Distributional equivalence and subcompositional coherence in the analysis of compositional data, contingency tables and ratio-scale measurements.  _J Classif_ , 26:29–54, 2009.
  * Greenacre et al. [2020] M. Greenacre, E. Grunsky, and J. Bacon-Shone.  A comparison of amalgamation and isometric logratios in compositional data analysis.  _Comput Geosc_ , 148:104621, 2020.
  * R Core Team [2021] 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/.
  * Ward [1963] J. H. Ward.  Hierarchical grouping to optimize an objective function.  _J Amer Stat Assoc_ , 58:236–244, 1963.


# Nonlinear Eigen-approach ADMM for Sparse Optimization on Stiefel Manifold

Jiawei Wang  Department of Mathematics, Hong Kong Baptist University  Rencang
Li  Department of Mathematics, University of Texas at Arlington  Richard Yi Da
Xu  Department of Mathematics, Hong Kong Baptist University

(May 28, 2023)

###### Abstract

With the growing interest and applications in machine learning and data
science, finding an efficient method to sparse analysis the high-dimensional
data and optimizing a dimension reduction model to extract lower dimensional
features has becoming more and more important. Orthogonal constraints (Stiefel
manifold) is a commonly met constraint in these applications, and the sparsity
is usually enforced through the element-wise L1 norm. Many applications can be
found on optimization over Stiefel manifold within the area of physics and
machine learning. In this paper, we propose a novel idea by tackling the
Stiefel manifold through an nonlinear eigen-approach by first using ADMM to
split the problem into smooth optimization over manifold and convex non-smooth
optimization, and then transforming the former into the form of nonlinear
eigenvalue problem with eigenvector dependency (NEPv) which is solved by self-
consistent field (SCF) iteration, and the latter can be found to have an
closed-form solution through proximal gradient. Compared with existing
methods, our proposed algorithm takes the advantage of specific structure of
the objective function, and has efficient convergence results under mild
assumptions.

Keywords: Stiefel manifold, Nonlinear eigenvalue problem with eigenvector
dependency (NEPv), self-consistent field (SCF) iteration, ADMM, proximal
algorithm, orthogonal constraints

##  1 Introduction

###  1.1 Problem Description

In this paper, we are concerned with the following non-smooth composite
minimization problem

$\begin{split}\min_{X}&\;\;f(X)+r(X)\\\
\mathrm{s.t.}&\;\;X\in\mathcal{S}_{n,p}\end{split}$ (1.1)

where $\mathcal{S}_{n,p}$ stands for the Stiefel manifold
$\mathcal{S}_{n,p}=\\{X\in\mathbb{R}^{n\times p}\mid X^{\top}X=I_{p}$},
$f:\mathbb{R}^{n\times p}\mapsto\mathbb{R}$ needs to satisfy the following
Assumption 1.1, and $r:\mathbb{R}^{n\times p}\mapsto\mathbb{R}$ represents the
convex but possibly non-smooth function whose proximal operator can be easily
evaluated.

######  Assumption 1.1.

$f(X)$ is differentiable and possibly non-convex, and $\nabla f(X)$ is
Lipschitz continuous with Lipschitz constant $L_{f}$.

In many applications, $r(\cdot)$ is usually being taken as the
$\ell_{1}$-related regularization term to enforce sparsity, and the frequently
encountered ones include for example the element-wise matrix $\ell_{1}$-norm
and the $\ell_{m,1}$-norm

$r(X)\coloneqq\begin{cases}\mu\|X\|_{1}=\sum_{ij}\mu|X_{ij}|\\\ \|\Gamma
X\|_{m,1}=\sum_{i=1}^{n}\gamma_{i}\|X_{i\cdot}\|_{m}\end{cases}$ (1.2)

where $\mu>0$, $\gamma_{i}\geq 0$ are weighting parameters and $\|\cdot\|_{m}$
denotes the $\ell_{m}$ vector norm for $X_{i\cdot}$, the $i$th row of matrix
$X$.

###  1.2 Existing Methods

Optimization on Stiefel manifold has been a hot research topic in the past few
years. [1] gives a very comprehensive and structured study on optimization
method over smooth manifold. However, most of the existing methods
[2][3][4][5][6] are designed for smooth objective function by using the
derivative of objective function and thus cannot be applied to the non-smooth
objective function such as problem (1.1).

Due to the non-convexity, non-smoothness as well as the manifold constraint,
algorithms that can be applied here is limited, and the existing ones can be
basically divided into two categories. The first type focuses on using
Riemannian algorithm to handle the manifold constraint [7][8][9][10].
Riemannian subgradient (1.3) is used in [7][8] to deal with non-smoothness,

$X_{k}={\rm Retr}_{X_{k-1}}(-\alpha_{k-1}\nu_{k-1})$ (1.3)

where $\alpha_{k-1}$ is the stepsize, $\nu_{k-1}$ denotes the Riemannian
subgradient, and Retr is the retraction mapping. However, the efficiency of
these proposed algorithm is not ideal as the subgradient method has relatively
slow convergence rate, and the method is too general where it does not utilize
the unique properties from the problem. The idea of MADMM [9] is relatively
close to our algorithm where the variable splitting and ADMM is applied first
to divide the original minimization into two sub-problems, but then an inner
algorithm of Riemannian gradient is used to solve the smooth manifold
constraint sub-problem. This splitting technique is not new and it is first
proposed in [11] by Lai and Osher where it is called splitting orthogonal
constraints (SOC). However, there are no convergence guarantees for both SOC
based method in [11] as well as MADMM. Chen et al. [10] proposed an
retraction-based proximal gradient called manifold proximal gradient method
(ManPG), where the proximal method is extended from the Euclidean setting to
manifold setting. In ManPG, after computing the proximal algorithm on tangent
space by an inner semi-smooth Newton algorithm, a retraction is then adopted
to map the iteration points back to the smooth manifold. Here, the
computational cost for solving the proximal mapping sub-problem is high, which
may not be efficient in certain problems.

The methods within the second category apply the Lagrangian based approach
(e.g. augmented Lagrangian multiplier method) or indicator function on the
orthogonal constraint to transform the minimization (1.1) unconditional. In
[12], Chen et al. proposed an proximal alternating minimized augmented
Lagrangian (PAMAL) method to solve some specific $\ell_{1}$-regularized
optimization problem such as the compressed modes [13][14] in physics or the
feature selection problem [15][16] with minor modifications. The PAMAL is
based on an proximal regularized three-block Gauss-Seidel method applied to
the augmented Lagrangian function 1.4 after denoting the orthogonal constraint
using an indicator function and splitting the variable twice.

$\begin{split}\mathcal{L}(X,Y,Z,\Lambda^{1},\Lambda^{2};\beta)=&{\rm
trace}(X^{T}HX)+\frac{1}{\mu}\|Y\|_{1}+\delta_{\mathcal{S}}(Z)+\langle\Lambda^{1},Y-X\rangle\\\
&+\frac{\beta}{2}\|Y-X\|_{F}^{2}+\langle\Lambda^{2},Z-X\rangle+\frac{\beta}{2}\|Z-X\|_{F}^{2}\end{split}$
(1.4)

where $X,Y,Z,\Lambda^{i}\in\mathbb{R}^{n\times p}$ for $i=1,2$, $\mu,\beta>0$
and $H$ is the (discrete) Hamiltonian. The update of the primal variable
$X_{k},Y_{k},Z_{k}$ is up to some predefined tolerance level
$\|\Theta_{k}\|_{\infty}\leq\epsilon_{k},\Theta_{k}\in\partial\mathcal{L}(X_{k},Y_{k},Z_{k},\Lambda^{1}_{k},\Lambda^{2}_{k};\beta_{k})$.
Apart from the relatively complicated algorithm structure, the result of PAMAL
is subjected and sensitive to its parameters settings like the three proximal
parameters. When $r(\cdot)$ within (1.1) is the $\ell_{m,1}$ norm ($m>1$) ,
[17] proposes an proximal gradient method for penalty models with a compact
and convex constraints (PenCPG) based on using augmented Lagrangian method on
orthogonal constraint and instead of dual ascent, the authors in [17] express
the Lagrange multiplier through an close-form expression of the primal
variable:

$\min_{X\in\mathcal{M}}g(X):=f(X)+r(X)-\frac{1}{2}\langle\Lambda(X),X^{T}X-I_{p}\rangle+\frac{\beta}{4}\|X^{T}X-I_{p}\|_{F}^{2}$
(1.5)

where $\mathcal{M}\subseteq\mathbb{R}^{n\times p}$ denotes a compact convex
set includes the Stiefel manifold. Although this method gets rid of the dual
variable and utilize the fast convergence rate of proximal algorithm, the
limitation is that the closed-form expression for $\Lambda$ only exist for
$r(\cdot)$ being $\ell_{m,1}$ norm ($m>1$).

###  1.3 Motivation and Contribution

Inspired by the novel idea of nonlinear eigenvalue problem with eigenvector
dependency proposed by the authors in [18], algorithms based on the NEPv
formulation and then SCF iteration are designed to solve the problems like
$\theta$-trace ratio [6] and orthogonal canonical correlation analysis (OCCA)
[19] which are defined on Stiefel manifold arise. However, the partial
derivative of the objective function and its corresponding gradient on
manifold are required during the process of transforming the first-order
optimality condition (KKT) into the NEPv structure, which means the problem
considered needs to be smooth. Therefore, ADMM is considered on problem (1.1)
to split the difficulty of the original minimization for which the smoothness
condition for the construction of NEPv can then be satisfied.

Many of the optimization problems in machine learning defined over Stiefel
manifold (e.g. principal component analysis (PCA), orthogonal linear
discriminant analysis (OLDA), canonical correlations analysis (CCA) etc.) are
aim at (linear) dimensionality reduction, which is the cornerstone when
studying the high dimensional data, and the related function can be expressed
through the matrix trace function. The $p$ in problem (1.1) is usually
referred to as the number of selected features and is a very small number
compared to $n$ ($p\ll n$), the original dimensionality from data points. The
NEPv approach from our NEPvADMM can be designed to better take advantage of
these circumstances as well as the structure of each specific objective
function to further lower the computational cost than the Riemannian gradient
descent method from MADMM.

###  1.4 Notation and Organization

In this paper, $\mathbb{R}^{n\times p}$ is the set of $n\times p$ real
matrices, $\mathbb{S}^{n\times n}$ denotes all $n\times n$ real symmetric
matrices. $I_{n}$ is the identity matrix and $X^{T}$ is the transpose of a
matrix. The Euclidean inner product of $X,Y\in\mathbb{R}^{n\times p}$ is
defined through the trace function $\langle X,Y\rangle={\rm tr}(X^{T}Y)$,
where ${\rm tr}(H)=\sum_{i}H_{ii}$ for $H\in\mathbb{R}^{n\times n}$. Here we
will use $\lambda_{i}(A)$ to denote the $i$-th largest eigenvalue of matrix
$A\in\mathbb{R}^{n\times n}$ and $\sigma_{i}(A)$ being its $i$-th largest
singular value.

$\|A\|_{F}=(\sum_{i,j}A^{2}_{i,j})^{\frac{1}{2}}=(\sum_{j}\sigma^{2}_{j})^{\frac{1}{2}},$
$\|A\|_{2}=\max_{j}\sigma_{j},$ $\|A\|_{\rm trace}=\sum_{j}\sigma_{j},$

denote the Frobenius norm, matrix spectral norm, and matrix trace norm
respectively. Finally, $A\in\mathbb{S}^{n\times n}\succeq(\succ)0$ means the
matrix is positive semi-definite (positive definite).

The rest of the paper is organized as follows. In Section 2, we will explain
the methodology of the NEPv ADMM method in detail and provide its algorithm.
And in Section 3, real-life applications will be introduced and relative
numerical experiments are implemented. Finally, we draw the conclusion in
Section 4.

##  2 NEPv ADMM

The objective function 1.1 is a combination of a smooth (possibly non-convex)
function and a convex non-smooth function, it is natural to think of using
alternating direction method of multipliers (ADMM) to split the difficulty of
the original problem and handle each sub-problem separately, where one is
smooth optimization on manifold and the other is simply convex optimization.

First, we split the variables and transform the original problem 1.1 into an
equivalent one

$\begin{split}\min_{X,Y}&\;\;-\phi(X)+r(Y)\\\
\mathrm{s.t.}\;\;&X\in\mathcal{S}_{n,p},Y\in\mathbb{R}^{n\times
p},X=Y\end{split}$ (2.1)

where $\phi(X)=-f(X)$.

The augmented Lagrangian function is of 2.1 is then

$\mathcal{L}_{\beta}(X,Y;\Lambda)=-\phi(X)+r(Y)+\langle\Lambda,X-Y\rangle+\frac{\beta}{2}\|X-Y\|_{F}^{2},$
(2.2)

Applying ADMM on the augmented Lagrangian function and update each variable
iteratively

$\begin{split}X_{k+1}&:=\mathop{\arg\max}\limits_{X\in\mathcal{M}}-\mathcal{L}_{\beta}(X,Y_{k};\Lambda_{k})\\\
Y_{k+1}&:=\mathop{\arg\min}\limits_{Y}\mathcal{L}_{\beta}(X_{k+1},Y;\Lambda_{k})\\\
\Lambda_{k+1}&:=\Lambda_{k}+\beta(X_{k+1}-Y_{k+1}).\end{split}$ (2.3)

From 2.3, due to the existence of the Stiefel manifold constraints, the X sub-
problem is an optimization on smooth manifold. Instead of using the commonly
applied methods such as Riemannian gradient descent, we apply a more
customized nonlinear eigenvalue problem with eigenvector dependency (NEPv)
approach to tackle it which is more efficient as $f(X)$ within 1.1 is usually
about dimensionality reduction. The general idea of the NEPv approach is to
get the first order optimality condition by taking the gradient of the
objection function first, and then transform the equality into the form of
nonlinear eigenvalue problem with eigenvector dependency and followed by using
the self-consistent field (SCF) iteration.

Denoting $g_{k}(X)=-\mathcal{L}_{\beta}(X,Y_{k};\Lambda_{k})$ for simplicity,
the gradient of $g_{k}(X)$ on Stiefel manifold is

${\rm grad}\;g_{k}(X)=\frac{\partial g_{k}(X)}{\partial X}-X{\rm
sym}(X^{T}\frac{\partial g_{k}(X)}{\partial X})$ (2.4)

where ${\rm sym}(A)=\frac{A+A^{T}}{2}$.

Therefore, the KKT condition of the X sub-problem can be transformed into

$\frac{\partial g_{k}(X)}{\partial X}=X\hat{\Omega}$ (2.5)

where $X\in\mathbb{O}^{n\times p}$ and
$\hat{\Omega}^{T}=\hat{\Omega}\in\mathbb{R}^{p\times p}$.

The partial derivative $\frac{\partial g_{k}(X)}{\partial X}$ is computed by
treating each entries within $X$ independent from each other:

$\frac{\partial g_{k}(X)}{\partial X}=\frac{\partial\phi(X)}{\partial
X}-\Lambda_{k}-\beta(X-Y_{k})$ (2.6)

Here we do not have the specific form of $\phi(X)$, thus we take the general
method to transform the left hand side of 2.5 into the form of NEPv $H(X)X$
where $H(X)\in\mathbb{R^{n\times n}}$ is a symmetric matrix

$\begin{split}H(X)&=\frac{\partial(g_{k}(X)+\beta X)}{\partial
X}X^{T}+X(\frac{\partial(g_{k}(X)+\beta X)}{\partial X})^{T}-\beta I_{n}\\\
&=\frac{\partial\phi(X)}{\partial X}X^{T}+X(\frac{\partial\phi(X)}{\partial
X})^{T}+D_{k}X^{T}+XD_{k}^{T}-\beta I_{n}\end{split}$ (2.7)

where $D_{k}=\beta Y_{k}-\Lambda_{k}$.

But special and simpler formulations are often available if the specific form
of $\phi(X)$ is known. For example, in sparse PCA, $\phi(X)=\frac{1}{2}{\rm
trace}(X^{T}AX)$ where $A\succeq 0$, we can take $H(X)=A-\beta
I_{n}+D_{k}X^{T}+XD_{k}$.

The NEPv of the X sub-problem takes the form

$H(X)X=X\Omega,\;X\in\mathcal{S}_{n,p},\;\Omega^{T}=\Omega\in\mathbb{R}^{p\times
p}$ (2.8)

######  Theorem 1.

$X\in\mathcal{S}_{n,p}$ satisfies the first order optimality condition of
$\mathop{\max}\limits_{X\in\mathcal{S}_{n,p}}g_{k}(X)$ if and only if it is a
solution of NEPv (2.8) and $X^{T}(D_{k}+\frac{\partial\phi(X)}{\partial X})$
is symmetric.

To make sure the monotonic increase of each iteration step with the X sub-
problem, we assume the NEPv Ansatz holds, meaning for
$\hat{X}\in\mathcal{S}_{n,p}$, $X\in\mathbb{X}\subseteq\mathcal{S}_{n,p}$, if

${\rm trace}(\hat{X}^{T}H(X)\hat{X})\geq{\rm trace}(X^{T}H(X)X)+\eta\;\;\;{\rm
for\;some}\;\eta\in\mathbb{R}$ (2.9)

Then there exist $Q\in\mathcal{S}_{p,p}$ with
$\tilde{X}=\hat{X}Q\in\mathbb{X}$ such that $g_{k}(\tilde{X})\geq
g_{k}(X)+\omega\eta$.

Next, the SCF iteration described in Algorithm 1 can be used

Algorithm 1 Self-Consistent Field Iteration

$X_{0}\in\mathbb{X}$, set $k=0$

while certain stopping criterion is not reached do Construct $H_{k}=H(X_{k})$;
Compute orthonormal eigenbasis matrix $\hat{X}_{k+1}$ associated with the $p$
largest eigenvalues of $H_{k}$; Compute $Q_{k+1}\in\mathbb{X}$ and set
$X_{k+1}=\hat{X}_{k+1}Q_{k+1}\in\mathbb{X}$; Set $k\coloneqq k+1$.

end while

Return the last $X_{k}$ as a maximizer of $g_{k}(X)$.

Since the exact optimizer of X sub-problem is hard to solve, we proposed to
replace it by only one SCF iteration in the X sub-problem.

The Y sub-problem is a convex optimization, and it has a closed-form solution
through proximal gradient method.

######  Lemma 1.

The Y sub-problem has closed-form solution $Y_{k+1}={\rm
Prox}_{\frac{r(\cdot)}{\beta}}(X_{k+1}+\frac{\Lambda_{k}}{\beta})$, and we
have $\Lambda_{k+1}\in\partial r(Y_{k+1})$. For $r(\cdot)$ being the element-
wise $\mathcal{L}_{1}$-norm, the proximal operator can be solved by soft-
thresholding function.

###### Proof.

The Y sub-problem:

$Y_{k+1}:=\mathop{\arg\min}\limits_{Y\in\mathbb{R}^{n\times
p}}r(Y)+\langle\Lambda_{k},X_{k+1}-Y\rangle+\frac{\beta}{2}\|X_{k+1}-Y\|_{F}^{2}$
(2.10)

It has optimality condition:

$\begin{split}0&\in\partial r(Y_{k+1})-\Lambda_{k}+\beta(Y_{k+1}-X_{k+1})\\\
0&\in\partial\frac{1}{\beta}r(Y_{k+1})+Y_{k+1}-(X_{k+1}+\frac{1}{\beta}\Lambda_{k})\end{split}$
(2.11)

which can be expressed through a one-step proximal gradient: $Y_{k+1}={\rm
Prox}_{\frac{r(\cdot)}{\beta}}(X_{k+1}+\frac{\Lambda_{k}}{\beta})$.

Rearranging the terms of the first equation within (2.11):

$\begin{split}\Lambda_{k}+\beta(X_{k+1}-Y_{k+1})&\in\partial r(Y_{k+1})\\\
\Lambda_{k+1}&\in\partial r(Y_{k+1})\end{split}$ (2.12)

∎

Our NEPv based ADMM can then be summarized by the following Algorithm 2, where
we can see the purcedure is straightforward and no inner loops are included.

Algorithm 2 NEPv based ADMM Algorithm

$X_{0}\in\mathbb{X}\subseteq\mathcal{S}_{n,p}$, $Y_{0}$, $\beta$, set $k=0$

while certain stopping criterion is not reached do Construct $H_{k}=H(X_{k})$;
Compute orthonormal eigenbasis matrix $\hat{X}_{k+1}$ associated with the $p$
largest eigenvalues of $H_{k}$; Compute $Q_{k+1}\in\mathbb{X}$ and update
$X_{k+1}=\hat{X}_{k+1}Q_{k+1}\in\mathbb{X}$; Update $Y_{k+1}:={\rm
Prox}_{\frac{r(\cdot)}{\beta}}(X_{k+1}+\frac{\Lambda_{k}}{\beta})$; Update
$\Lambda_{k+1}:=\Lambda_{k}+\beta(X_{k+1}-Y_{k+1})$; Set $k\coloneqq k+1$.

end while

For standard two-block ADMM algorithm which has already been proven to
converge to the first-order optimality point, the two functions $f(X)$ and
$r(X)$ from the composite objective function both need to be proper, closed
convex functions, as well as each of the two sub-problems are optimized
exactly. In our ADMM structure from Algorithm 2, the X sub-problem is non-
convex and is an inexact optimization. Although several works have expanded
the convergence condition of ADMM to some mild assumption on the objective
function, none of them cover the type of the problem considered here.

##  3 Numerical Experiments

###  3.1 Applications

Many applications can be found suited to the problem formulation (1.1) within
the area of machine learning, data science, and physics. We list some of them
in the following.

Example 4.1 (sparse principle component analysis (sPCA)) Principal component
analysis is one of the most widely used unsupervised learning method in
machine learning for dimension reduction. It was first proposed by Pearson in
1901 [20] as a method to minimize the projection residuals, but nowadays it is
often treated as finding a projection that can maximize the variance, and the
two perspectives are equivalent. However, in traditional PCA, it is difficult
to interpret the corresponding principal components as the each of them is a
dense linear combination of all the original features [21]. Nowadays, massive
amounts of data needs to be processed and these data are of high-dimensional
(usually the number of features is larger than the number of data
observations) in general, and the lack of feature selection has restricted the
performance of traditional PCA in many applications such as the sparse coding
of image processing [22][23] as well as the natural language processing [24].

To address these issues by improving the interpretability as well as reducing
the storage required [25], sparse PCA was proposed using the idea that most of
the variation in the data can be explained using only a small number of
significant original features and thus the principal components are formed by
sparse combinations instead. There are different ways to enforce sparsity,
here we consider adding the $\ell_{1}$-regularization term to the original PCA
objective function as it is generally easy to solve and can be useful in many
applications [26]. The problem formulation is as follows (3.1).

$\begin{split}\min_{X}&-\frac{1}{2}tr(X^{\top}AA^{\top}X)+\mu\|X\|_{1}\\\
&\mathrm{s.t.}\;\;\;X\in\mathcal{S}_{n,p}\end{split}$ (3.1)

where $A$ is $\mathbb{R}^{n\times m}$ data matrix and $AA^{T}\succeq 0$. $n$
denote the dimension of a single data, meaning the types of features, $m$ is
the amount of data observations, and $\mu>0$ here is a weighting parameter for
the problem and (3.1) reduces back to PCA if $\mu=0$ (the leading $p$
eigenvalues of $AA^{T}$).

Example 4.2 (orthogonal dictionary learning (ODL)) Given the data matrix
$A\in\mathbb{n\times p}$ with the column vectors represent the data points
$\mathbf{a}_{i},i=1,...,p$ and $p>>n$, the purpose of ODL is to find an
orthogonal matrix $X\in\mathbb{R}^{n\times n}$ so that $A$ can be compactly
represented [27][28][29]. The problem has the following generalized
mathematical format:

$\begin{split}\min_{X}&\;\;\|A^{T}X\|_{1}\\\ {\rm
s.t.}&\;\;X\in\mathcal{S}_{n,n}\end{split}$ (3.2)

The $A^{T}X$ from the objective function of (3.2) can be seen as the
correlation matrix, which we want it to be sparse by controlling the
orthogonal basis matrix $X$, and the element-wise $\ell_{1}$ norm here is used
to enforce the sparsity. Based on (3.2), [30] introduces an application called
robust subspace recovery (RSR), and the corresponding dual principal component
pursuit (DPCP) approach for solving RSR can be found in [31][32][33], where
they transform the problem into:

$\begin{split}\min_{X}&\;\;\|A^{T}X\|_{1}\\\
s.t.&\;\;X\in\mathcal{S}_{n,p}\end{split}$ (3.3)

where $p=n-d$ and $d$ is the dimensionality of the subspace spanned by the
inlier point from the data matrix $A$.

Example 4.3 (compressed modes in physics) For the independent-particle
Schrödinger equation about a finite system of electrons, the compressed modes
problem is to compute a spatially localized solutions to its eigenvalue
problem, where the sparsity is enforced through adding the $\ell_{1}$
regularization term after the wave function. According to [34], the problem
has the following formulation:

$\begin{split}\min_{X\in\mathbb{R}^{n\times p}}&\;\;{\rm
trace}(X^{T}HX)+\frac{1}{\mu}\|X\|_{1}\\\ {\rm
s.t.}&\;\;X^{T}X=I_{p}\end{split}$ (3.4)

where $H\in\mathbb{R}^{n\times n}$ represents the (discrete) Hamiltonian,
parameter $\mu>0$ is predefined for balancing the weight of sparsity.

Example 4.4 (sparse spectral clustering) Spectral clustering (SC) is one of
the most commonly used data clustering technique, and it requires to find out
a lower dimensional embedding $X$ of the original data through calculating the
eigenvalues of the normalized Laplacian matrix first, then the final results
can be derived through applying the K-means algorithm on $X^{T}$. The sparse
SC is an extension of SC by implementing the idea of sparse regularization
into SC, and according to the authors from [35], the mathematical expression
of the problem can be formulated as following:

$\begin{split}\min_{Y\in\mathbb{R}^{n\times n},X\in\mathbb{R}^{n\times
p}}&\;\;\langle L,XX^{T}\rangle+g(Y)\\\ {\rm
s.t.}&\;\;Y=XX^{T},X^{T}X=I.\end{split}$ (3.5)

where $L$ denotes the normalized Laplacian matrix and $g:\mathbb{R}^{n\times
n}\rightarrow\mathbb{R}$ is the sparse regularizer which may be selected
differently depending on problem settings.

Example 4.5 (unsupervised feature selection) Feature selection (FS) is about
getting rid of the noisy features and selecting a feature subset with much
smaller dimensionality from some high-dimensional feature set, so that the
researcher may have a better data representation. In contrast to supervised FS
such as trace ratio and Fisher score, unsupervised FS is much more difficult
to handle because of its lack of label information. By adopting the assumption
that a linear classifier can be used to predict the class label of input data,
[16] incorporates the $\ell_{2,1}$ norm and discriminative analysis to model
the unsupervised FS as a manifold optimization problem in the form of problem
(1.1):

$\begin{split}\min_{W}&\;\;{\rm trace}(W^{T}MW)+\gamma\|W\|_{2,1}\\\ {\rm
s.t.}&\;\;W\in\mathcal{S}_{n,p}\end{split}$ (3.6)

where $M\in\mathbb{S}^{n\times n}$ is calculated from the input data.

###  3.2 Numerical Experiments on sparse PCA

In this section, the efficiency of our NEPvADMM will be illustrated and
presented through conducting numerical experiments on the problem of sparse
PCA. Firstly, the test problem and the basic setting such as the stopping
criterion and default parameter values will be explained. Next, we will
compare the numerical performances between our algorithm with the existing
ones and draw some conclusions. All the algorithms are performed in MATLAB
R2022b, and the experiments are conducted on a laptop with Intel Core
i5-1035G1 CPU @ 1.00GHz and 8G of RAM.

####  3.2.1 Experimental Settings

To implement the NEPvADMM, the variables within the original sparse PCA (3.1)
are splitted and transformed to

$\begin{split}\min_{X,Y}&\;\;-\frac{1}{2}{\rm trace}(X^{T}MX)+\mu\|Y\|_{1}\\\
{\rm s.t.}&\;\;X=Y,X\in\mathcal{S}_{n,p}\end{split}$ (3.7)

where $M=A^{T}A\in\mathbb{R}^{n\times n}$. The augmented Lagrangian function
of (3.7) is

$\mathcal{L}_{\beta}(X,Y;\Lambda)=-\frac{1}{2}{\rm
trace}(X^{T}MX)+\mu\|Y\|_{1}+\langle\Lambda,X-Y\rangle+\frac{\beta}{2}\|X-Y\|_{F}^{2},$
(3.8)

After applying the ADMM (2.3) to (3.8), we need to check the assumption that
the function value (3.9) from $X$-subproblem is monotonic increasing under the
SCF iteration

$g_{k}(X)=\frac{1}{2}{\rm
trace}(X^{T}MX)-\langle\Lambda_{k},X-Y_{k}\rangle-\frac{\beta}{2}\|X-Y_{k}\|_{F}^{2}$
(3.9)

The following Lemma tells us whether the NEPv ansatz holds for sparse PCA:

######  Lemma 2.

Given $X\in\mathcal{S}_{n,p}$, the NEPv ansatz holds for the objective
function (3.9) of X-subproblem from sparse PCA.

###### Proof.

First, we need to transform the gradient of (3.9) into the form of NEPv to
contract the matrix $H(X)$. In Section 2, for the sparse PCA, a simpler form
of $H(X)$ can be derived instead of directly applying (2.7)

$H(X)=M-\beta I_{n}+D_{k}X^{T}+XD_{k}$ (3.10)

where $D_{k}=\beta Y_{k}-\Lambda_{k}$. Next, for
$X,\hat{X}\in\mathcal{S}_{n,p}$, if

${\rm trace}(\hat{X}^{T}H(X)\hat{X})\geq{\rm trace}(X^{T}H(X)X),$ (3.11)

Inserting (3.10) into (3.11), through simple derivation, the left-hand side of
(3.11) is

$\begin{split}{\rm trace}(\hat{X}^{T}H(X)\hat{X})&={\rm
trace}(\hat{X}^{T}M\hat{X})-\beta{\rm trace}(\hat{X}^{T}\hat{X})+2{\rm
trace}(\hat{X}^{T}D_{k}X^{T}\hat{X})\\\ &\leq{\rm
trace}(\hat{X}^{T}M\hat{X})-\beta{\rm
trace}(\hat{X}^{T}\hat{X})+\|\hat{X}^{T}D_{k}\|_{\rm trace}\\\ &={\rm
trace}(\tilde{X}^{T}M\tilde{X})-\beta{\rm trace}(\tilde{X}^{T}\tilde{X})+{\rm
trace}(\tilde{X}^{T}D_{k})\end{split}$ (3.12)

where the inequality within (3.12) is done by von Neumann’s trace inequality
and $\tilde{X}=\hat{X}(UV^{T})$ in the last equality is defined through
singular value decomposition (SVD) $\hat{X}^{T}D_{k}=U\Sigma V^{T}$.
Similarly, the right-hand side of (3.11) has the form of

${\rm trace}(X^{T}H(X)X)={\rm trace}(X^{T}MX)-\beta{\rm trace}(X^{T}X)+2{\rm
trace}(X^{T}D_{k})$ (3.13)

Adding the term ${\rm trace}(Y_{k}^{T}\Lambda_{k})-\beta{\rm
trace}(Y_{k}^{T}Y_{k})$ on both side of (3.12) and (3.13), the inequality
(3.11) can be transformed into the form of (3.9)

$g_{k}(\tilde{X})\geq g_{k}(X)$ (3.14)

which shows that the ansatz is satisfied. ∎

The data matrix $A\in\mathbb{R}^{n\times p}:=\mathbf{randn}(n,p)$ is randomly
generated, where $n$ is chosen from $\\{500,1000\\}$, $p$ from $\\{5,10\\}$.
We set $\beta=20$ and $\mu$ from $\\{0.5,1\\}$. The stopping criterion is
defined as

$\frac{|F(Y_{k+1})-F(Y_{k})|}{|F(Y_{k})|}<10^{-8}$

where $F(\cdot):=f(\cdot)+r(\cdot)$.

####  3.2.2 Numerical Results

We compare the performance of our NEPvADMM with the Riemannian subgradient
descent method and the MADMM approach. The following Table 1 presents the
results from three aspects: the objection value and the sparisty, where the
sparsity formula is calculated by the percentage of the zero element within
the matrix. From the table, we can see that the Riemannian subgradient method
is unable to iteration a sparse result while the NEPvADMM and MADMM can, and
the sparisty are very high both the latter two methods. Besides, the objective
value generated by our approach and MADMM are very close, however, the ones
derived by RSG is high compared to the others.

Table 1: Results comparison between NEPvADMM, Riemannian subgradient method and MADMM from the aspects of objection value and sparsity Setting | NEPvADMM | RSG | MADMM
---|---|---|---
$\mu$ | (n,p) | Obj. | Spa. | Obj. | Spa. | Obj. | Spa.
$0.5$ | (300,5) | $1.1580$ | $0.9967$ | $3.0364$ | $0$ | $1.1786$ | $0.9967$
$0.5$ | (300,50) | $12.5907$ | $0.9967$ | $30.7831$ | $0$ | $12.5042$ | $0.9960$
$0.5$ | (500,5) | $1.2330$ | $0.9980$ | $4.3290$ | $0$ | $1.2330$ | $0.9980$
$0.5$ | (500,50) | $12.6830$ | $0.9980$ | $43.6577$ | $0$ | $12.6816$ | $0.9980$
$1$ | (300,5) | $3.7739$ | $0.9967$ | $10.7888$ | $0$ | $3.7739$ | $0.9967$
$1$ | (300,50) | $36.9932$ | $0.9967$ | $110.0859$ | $0$ | $37.1584$ | $0.9967$
$1$ | (500,5) | $3.7146$ | $0.9980$ | $11.5442$ | $0$ | $3.7146$ | $0.9980$
$1$ | (500,50) | $38.0765$ | $0.9980$ | $157.2606$ | $0$ | $38.1054$ | $0.9980$

To have a better presentation of the efficiency of these three methods, Figure
1 shows how the objective value changes along with iteration counts as well as
the CPU times. From Figure 1, the objective value of NEPvADMM and MADMM are
similar while the RSG (blue line) performs badly, and it is clear that our
approach requires much lower iteration counts as well as CPU time to converge
in both four cases.

(a) n = 500, p = 5

(b) n = 500, p = 50

(c) n = 500, p = 5

(d) n = 500, p = 50

Figure 1: Comparison of the iteration counts required and the CPU time (in
seconds) consumed among the NEPvADMM, Riemannian gradient method, and the
MADMM approach for $\mu=1$

##  4 Conclusion

In this paper, we introduced an innovative method to address the non-smooth
composite optimization on Stiefel manifold by employing a nonlinear eigen-
approach (NEPv) together with alternating direction method of multipliers,
which fully utilize the unique properties from the problem such as dimension
reduction as well as sparsity, has great advantages compared to the limited
existing algorithms, and can be applied to many applications from machine
learning and data science. The steps from our algorithm in straightforward and
easy to implement. It can be customized to specific problems to further
improve the efficiency, and requires much smaller iteration counts to converge
compared to other existing approaches. We have strong belief in our method for
solving the non-smooth problems with orthogonal constraint. In future work, we
will not only conducted detailed mathematical analysis on the convergence
property of the algorithm, but also study and generalize the conditions for
objective functions where the NEPv ansatz is guaranteed to hold.

## References

  * [1] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds.  Princeton University Press, 2008.
  * [2] N. Xiao, X. Liu, and Y.-x. Yuan, “A class of smooth exact penalty function methods for optimization problems with orthogonality constraints,” Optimization Methods and Software, vol. 37, no. 4, pp. 1205–1241, 2022.
  * [3] B. Gao, X. Liu, and Y.-x. Yuan, “Parallelizable algorithms for optimization problems with orthogonality constraints,” SIAM Journal on Scientific Computing, vol. 41, no. 3, pp. A1949–A1983, 2019.
  * [4] B. Gao, X. Liu, X. Chen, and Y.-x. Yuan, “A new first-order algorithmic framework for optimization problems with orthogonality constraints,” SIAM Journal on Optimization, vol. 28, no. 1, pp. 302–332, 2018.
  * [5] R.-C. Li, “A theory of the nepv approach for optimization on the stiefel manifold,” arXiv preprint arXiv:2305.00091, 2023.
  * [6] L. Wang, L.-H. Zhang, and R.-C. Li, “Trace ratio optimization with an application to multi-view learning,” Mathematical Programming, vol. 201, no. 1, pp. 97–131, 2023.
  * [7] O. Ferreira and P. Oliveira, “Subgradient algorithm on riemannian manifolds,” Journal of Optimization Theory and Applications, vol. 97, pp. 93–104, 1998.
  * [8] P. Grohs and S. Hosseini, “$\varepsilon$-subgradient algorithms for locally lipschitz functions on riemannian manifolds,” Advances in Computational Mathematics, vol. 42, pp. 333–360, 2016.
  * [9] A. Kovnatsky, K. Glashoff, and M. M. Bronstein, “Madmm: a generic algorithm for non-smooth optimization on manifolds,” in Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part V 14, pp. 680–696, Springer, 2016.
  * [10] S. Chen, S. Ma, A. Man-Cho So, and T. Zhang, “Proximal gradient method for nonsmooth optimization over the stiefel manifold,” SIAM Journal on Optimization, vol. 30, no. 1, pp. 210–239, 2020.
  * [11] R. Lai and S. Osher, “A splitting method for orthogonality constrained problems,” Journal of Scientific Computing, vol. 58, pp. 431–449, 2014.
  * [12] W. Chen, H. Ji, and Y. You, “An augmented lagrangian method for 1-regularized optimization problems with orthogonality constraints,” SIAM Journal on Scientific Computing, vol. 38, no. 4, pp. B570–B592, 2016.
  * [13] V. Ozoliņš, R. Lai, R. Caflisch, and S. Osher, “Compressed modes for variational problems in mathematics and physics,” Proceedings of the National Academy of Sciences, vol. 110, no. 46, pp. 18368–18373, 2013.
  * [14] V. Ozoliņš, R. Lai, R. Caflisch, and S. Osher, “Compressed plane waves yield a compactly supported multiresolution basis for the laplace operator,” Proceedings of the National Academy of Sciences, vol. 111, no. 5, pp. 1691–1696, 2014.
  * [15] J. Tang and H. Liu, “Unsupervised feature selection for linked social media data,” in Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 904–912, 2012.
  * [16] Y. Yang, H. T. Shen, Z. Ma, Z. Huang, and X. Zhou, “$\ell$ 2, 1-norm regularized discriminative feature selection for unsupervised learning,” in IJCAI international joint conference on artificial intelligence, 2011.
  * [17] N. Xiao, X. Liu, and Y.-x. Yuan, “Exact penalty function for 2,1 norm minimization over the stiefel manifold,” SIAM Journal on Optimization, vol. 31, no. 4, pp. 3097–3126, 2021.
  * [18] Y. Cai, L.-H. Zhang, Z. Bai, and R.-C. Li, “On an eigenvector-dependent nonlinear eigenvalue problem,” SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 3, pp. 1360–1382, 2018.
  * [19] L.-H. Zhang, L. Wang, Z. Bai, and R.-C. Li, “A self-consistent-field iteration for orthogonal canonical correlation analysis,” IEEE transactions on pattern analysis and machine intelligence, vol. 44, no. 2, pp. 890–904, 2020.
  * [20] K. Pearson, “Liii. on lines and planes of closest fit to systems of points in space,” The London, Edinburgh, and Dublin philosophical magazine and journal of science, vol. 2, no. 11, pp. 559–572, 1901.
  * [21] H. Zou, T. Hastie, and R. Tibshirani, “Sparse principal component analysis,” Journal of computational and graphical statistics, vol. 15, no. 2, pp. 265–286, 2006.
  * [22] F. Chen and K. Rohe, “A new basis for sparse principal component analysis,” Journal of Computational and Graphical Statistics, pp. 1–14, 2023.
  * [23] A.-K. Seghouane, N. Shokouhi, and I. Koch, “Sparse principal component analysis with preserved sparsity pattern,” IEEE Transactions on Image Processing, vol. 28, no. 7, pp. 3274–3285, 2019.
  * [24] R. Drikvandi and O. Lawal, “Sparse principal component analysis for natural language processing,” Annals of data science, vol. 10, no. 1, pp. 25–41, 2023.
  * [25] L. Wang, X. Liu, and Y. Zhang, “A communication-efficient and privacy-aware distributed algorithm for sparse pca,” Computational Optimization and Applications, vol. 85, no. 3, pp. 1033–1072, 2023.
  * [26] N. Parikh, S. Boyd, et al., “Proximal algorithms,” Foundations and trends® in Optimization, vol. 1, no. 3, pp. 127–239, 2014.
  * [27] D. A. Spielman, H. Wang, and J. Wright, “Exact recovery of sparsely-used dictionaries,” in Conference on Learning Theory, pp. 37–1, JMLR Workshop and Conference Proceedings, 2012.
  * [28] Q. Qu, J. Sun, and J. Wright, “Finding a sparse vector in a subspace: Linear sparsity using alternating directions,” Advances in Neural Information Processing Systems, vol. 27, 2014.
  * [29] L. Demanet and P. Hand, “Scaling law for recovering the sparsest element in a subspace,” Information and Inference: A Journal of the IMA, vol. 3, no. 4, pp. 295–309, 2014.
  * [30] G. Lerman, M. B. McCoy, J. A. Tropp, and T. Zhang, “Robust computation of linear models by convex relaxation,” Foundations of Computational Mathematics, vol. 15, pp. 363–410, 2015.
  * [31] G. Lerman and T. Maunu, “An overview of robust subspace recovery,” Proceedings of the IEEE, vol. 106, no. 8, pp. 1380–1410, 2018.
  * [32] M. C. Tsakiris and R. Vidal, “Dual principal component pursuit,” Journal of Machine Learning Research, vol. 19, no. 18, pp. 1–50, 2018.
  * [33] Z. Zhu, Y. Wang, D. Robinson, D. Naiman, R. Vidal, and M. Tsakiris, “Dual principal component pursuit: Improved analysis and efficient algorithms,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [34] V. Ozoliņš, R. Lai, R. Caflisch, and S. Osher, “Compressed modes for variational problems in mathematics and physics,” Proceedings of the National Academy of Sciences, vol. 110, no. 46, pp. 18368–18373, 2013.
  * [35] C. Lu, J. Feng, Z. Lin, and S. Yan, “Nonconvex sparse spectral clustering by alternating direction method of multipliers and its convergence analysis,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 32, 2018.