# Symmetric Stair Preconditioning of Linear Systems
for Parallel Trajectory Optimization

Xueyi Bu1, Brian Plancher2 This material is based upon work supported by the
National Science Foundation (under Award 2246022). Any opinions, findings,
conclusions, or recommendations expressed in this material are those of the
authors and do not necessarily reflect those of the funding
organizations.1Xueyi Bu is with the School of Engineering and Applied Science,
Columbia University, New York, NY<EMAIL_ADDRESS>Plancher is with
Barnard College, Columbia University, New York, NY<EMAIL_ADDRESS>
###### Abstract

There has been a growing interest in parallel strategies for solving
trajectory optimization problems. One key step in many algorithmic approaches
to trajectory optimization is the solution of moderately-large and sparse
linear systems. Iterative methods are particularly well-suited for parallel
solves of such systems. However, fast and stable convergence of iterative
methods is reliant on the application of a high-quality preconditioner that
reduces the spread and increase the clustering of the eigenvalues of the
target matrix. To improve the performance of these approaches, we present a
new parallel-friendly symmetric stair preconditioner. We prove that our
preconditioner has advantageous theoretical properties when used in
conjunction with iterative methods for trajectory optimization such as a more
clustered eigenvalue spectrum. Numerical experiments with typical trajectory
optimization problems reveal that as compared to the best alternative parallel
preconditioner from the literature, our symmetric stair preconditioner
provides up to a 34% reduction in condition number and up to a 25% reduction
in the number of resulting linear system solver iterations.

##  I Introduction

Trajectory optimization algorithms [1] are a powerful, and state-of-the-art
set of tools for synthesizing dynamic motions for complex robots [2, 3, 4, 5,
6, 7]. There has been historical interest in parallel strategies [8] for
solving trajectory optimization problems and several more recent efforts have
shown that significant computational benefits are possible by exploiting
parallelism in different stages of the algorithm on heterogeneous hardware
platforms including multi-core CPUs [9, 10, 11, 12], GPUs [13, 14, 15, 16, 17,
18], and FPGAs [19, 20, 21, 22]. This shift toward parallelism is growing
increasingly important with the impending end of Moore’s Law and the end of
Dennard Scaling, which have led to a utilization wall that limits the
performance a single CPU chip can deliver [23, 24].

One key step in many algorithmic approaches to trajectory optimization is the
solution of moderately-large and sparse linear systems. These systems can be
solved by direct factorization or through iterative fixed point approaches.
While many state-of-the-art solvers leverage factorization-based approaches
[25, 26], iterative methods, like the Preconditioned Conjugate Gradient (PCG)
algorithm [27], are particularly well-suited for parallelism, as they are
computationally dominated by matrix-vector products and vector reductions [28,
29]. Furthermore, in the context of very-large-scale linear systems GPU
implementations of PCG, and other iterative methods, have already been shown
to outperform both factorization and iterative approaches on the CPU [30, 31].
However, fast and stable convergence of such iterative methods is reliant on
the use of a high-quality symmetric preconditioner that reduces the spread and
increase the clustering of the eigenvalues of the target matrix [28, 32, 33,
34].

To improve the performance of these approaches, we present a new parallel-
friendly symmetric stair preconditioner. We prove that our preconditioner has
advantageous theoretical properties when used in conjunction with iterative
methods for trajectory optimization such as a more clustered eigenvalue
spectrum than previous parallel preconditioners. Numerical experiments with
typical trajectory optimization problems reveal that as compared to the best
alternative parallel preconditioner from the literature, our symmetric stair
preconditioner provides up to a 34% reduction in condition number and up to a
25% reduction in the number of PCG iterations needed for convergence.

##  II Background and Related Work

###  II-A Direct Trajectory Optimization

Trajectory optimization [1], also known as numerical optimal control, solves
an (often) nonlinear optimization problem to compute a robot’s optimal path
through an environment as a series of states (x $\in\mathbb{R}^{n}$) and
controls (u $\in\mathbb{R}^{m}$). These problems assume a discrete-time
dynamical system,

$x_{k+1}=f(x_{k},u_{k},h),\quad x_{0}=x_{s},$ (1)

with a timestep $h$, and minimize an additive cost function,

$J(X,U)=\ell_{f}(x_{N})+\sum^{N-1}_{k=0}\ell(x_{k},u_{k}).$ (2)

One approach for solving such problems is through the use of _direct methods_.
At each iteration of a direct method, a quadratic approximation of the
nonlinear problem defined by Equations 1 and 2 is formed around a nominal
trajectory ($x,u$), resulting in a quadratic program (QP). The resulting KKT
system [35, 36] is then solved to update the nominal trajectory. This process
is repeated until convergence.

As the KKT system is a saddle point system, one approach to solving it is
through two step process using the symmetric positive definite111While $S$ is
actually symmetric negative definite, one can instead negate it and solve for
$-\lambda^{*}$, enabling us to simplify our proofs by treating it as symmetric
positive definite in the remainder of this work. _Schur Complement_ , $S$.
This process first solves for the optimal Lagrange multipliers, $\lambda^{*}$,
and then recovers the state and control update, $\delta z^{*}=[\delta
x^{*},\delta u^{*}]$ as follows (where $c$ and $C$ are the constraint value
and Jacobian, and $g$ and $G$ are the cost Jacobian and Hessian respectively):

$\begin{split}S&=-CG^{-1}C^{T}\hskip 20.0pt\gamma=c-CG^{-1}g\\\
S\lambda^{*}&=\gamma\hskip 54.0pt\delta z^{*}=G^{-1}(g-C^{T}\lambda^{*})\\\
\end{split}$ (3)

If we denote the dynamics Jacobians as $A_{k}=f_{x_{k}}$ and
$B_{k}=f_{u_{k}}$, and cost Jacobians and Hessians as $q_{k}=J_{x_{k}}$,
$r_{k}=J_{u_{k}}$, $Q_{k}=J_{x_{k}x_{k}}$, and $R_{k}=J_{u_{k}u_{k}}$, we can
define:

$\displaystyle\theta_{k}$
$\displaystyle=-A_{k}Q_{k}^{-1}A_{k}^{T}-B_{k}R_{k}^{-1}B_{k}^{T}-Q_{k+1}^{-1}$
(4) $\displaystyle\phi_{k}$ $\displaystyle=A_{k}Q_{k}^{-1}$
$\displaystyle\zeta_{k}$
$\displaystyle=A_{k}Q_{k}^{-1}q_{k}+B_{k}R_{k}^{-1}r_{k}-Q_{k+1}^{-1}q_{k+1}.$

The Schur complement then takes the following block-tridiagonal form:

$\displaystyle S=\left(\begin{matrix}-Q_{0}^{-1}&\phi_{0}^{T}&&&\\\
\phi_{0}&\theta_{0}&\phi_{1}^{T}&&\\\ &&\ddots&&\\\
&&\phi_{N-2}&\theta_{N-2}&\phi_{N-1}^{T}\\\
&&&\phi_{N-1}&\theta_{N-1}\end{matrix}\right)$ (5)
$\displaystyle\gamma=c+\left(-Q_{0}^{-1}q_{0}\quad\zeta_{0}\quad\zeta_{1}\quad\ldots\quad\zeta_{N-1}\right)^{T}.$

For most trajectory optimization problems $Q$ and $R$ are chosen to be
symmetric positive definite, or made symmetric positive definite by
construction (e.g., in the case of quadratic penalty methods), and $C^{T}$ has
full column rank. Furthermore, by construction $\theta$ is both invertible and
symmetric. Similarly, although $A_{k}$ and $B_{k}$ will vary depending on the
integrator in use, for most popular explicit and semi-implicit integrators
they take the form

$A=I+hM\quad B=hN,$ (6)

where $M$ and $N$ are functions of the system dynamics $f$. Additional
conditions, such as Lipschitz continuity, are usually imposed on $f$ to ensure
the local existence and uniqueness of the solution to the dynamical system
[37]. As a result, $\|M\|$ is bounded; thus, for sufficiently small $h$,
$h\|M\|<1$ and $A=I+hM$, and thus $\phi$, is invertible. Finally, given all of
the above, the dominant computational step in this approach is solving for
$\lambda^{*}$ in Equation 3.

###  II-B Parallel Iterative Solvers

There has been a significant amount of prior work developing general purpose
parallel iterative solvers, many of which target the GPU [28, 38, 39, 40, 41,
42, 43]. Iterative methods solve the problem $S\lambda^{*}=\gamma$ for a given
$S$ and $\gamma$ by iterative refining an estimate for $\lambda$ up to some
tolerance $\epsilon$. The most popular of these methods is the Conjugate
Gradient (CG) method, which is applicable to systems where $S$ is positive
definite. The convergence rate of CG is directly related to the spread of the
Eigenvalues of $S\in\mathbb{R}^{n\times n}$, converging faster when they are
clustered. We also note that clustered, and moderate in magnitude, eigenvalues
also have the added benefit of avoiding round-off and overflow errors caused
by iterative floating point math.

To improve the performance of CG, a symmetric preconditioning matrix
$\Phi\approx S$ is often applied to instead solve the equivalent problem
$\Phi^{-1}S\lambda^{*}=\Phi^{-1}\gamma$ [44]. The resulting preconditioned CG
(PCG) algorithm leverages matrix-vector products with $S$ and $\Phi^{-1}$, as
well as vector reductions, both parallel friendly operations. Finally, we note
that PCG requires a symmetric preconditioner [44].

###  II-C Parallel Preconditioners

There is a large literature on preconditioning [28, 45, 46, 47, 48, 49] and
there are many different preconditioners that are optimized for computation on
vector or parallel processors. The most popular of these are the Jacobi and
Block-Jacobi preconditioners. These set:

$\Phi=\text{diag}(S)\hskip 20.0pt\text{or}\hskip 20.0pt\Phi=\text{block-
diag}(S).$ (7)

Previous GPU based iterative solvers mainly leveraged these preconditioners
[42, 43]. For block banded matrices, alternating and overlapping block
preconditioners have also been used in previous work. These methods compute
$\Phi^{-1}$ as a sum of the inverse of block-diagonal matrices that compose
$S$ [28, 50]. Finally, Polynomial splitting preconditioners have also found
widespread usage [28, 34]. These follow the pattern $S=\Psi-P$ and compute a
preconditioner where:

$S^{-1}\approx\Phi^{-1}=(I+\Psi^{-1}P+(\Psi^{-1}P)^{2}\dots)\Psi^{-1}.$ (8)

We note that these preconditioners are only valid where:

$\rho(\Psi^{-1}P)<1,$ (9)

which is called a convergent splitting and will guarantee convergence when
used with the CG algorithm [28, 32, 33, 34]. Also, increasing the degree of
the polynomial computes a better approximation of $S$ and improves the
convergence rate of the resulting PCG algorithm. However, this requires more
computation to compute the preconditioner and also often creates a
preconditioner with a larger bandwidth, requiring more memory. For block-
banded matrices, like $S$ in our problem, the values in the true inverse decay
exponentially as one moves away from the diagonal [51], this creates a
tradeoff between the accuracy and both the memory and computational complexity
of the preconditioner.

##  III Stair Preconditioners

One polynomial splitting for (symmetric) block tridiagonal matrices that has
been shown to outperform Jacobi and Block-Jacobi methods is the stair based
splitting [52, 53]. This splitting comes in two types, left (type 1) and right
(type 2), depending upon the direction of the stair. For a 3x3 symmetric block
tridiagonal $S$:

$\begin{split}S=&\begin{bmatrix}D_{1}&O_{1}&0\\\
O_{1}^{T}&\;\;D_{2}\;\;&O_{2}\\\ 0&O_{2}^{T}&D_{3}\\\
\end{bmatrix}\end{split}$ (10)

The stair splittings are:

$\begin{split}\Psi_{l}=&\begin{bmatrix}{\color[rgb]{0,0,1}D_{1}}&{\color[rgb]{0,0.6,0}0}&0\\\
{\color[rgb]{0,0,1}O_{1}^{T}}&\;\;{\color[rgb]{0,0,1}D_{2}}\;\;&{\color[rgb]{0,0,1}O_{2}}\\\
0&{\color[rgb]{0,0.6,0}0}&{\color[rgb]{0,0,1}D_{3}}\\\ \end{bmatrix}\quad
P_{l}=-\begin{bmatrix}0&{\color[rgb]{0,0.6,0}O_{1}}&0\\\ \;\;0\;\;&0&\;0\;\\\
0&{\color[rgb]{0,0.6,0}O_{2}^{T}}&0\\\ \end{bmatrix}\\\
\Psi_{r}=&\begin{bmatrix}{\color[rgb]{0,0,1}D_{1}}&{\color[rgb]{0,0,1}O_{1}}&0\\\
{\color[rgb]{0,0.6,0}0}&\;\;{\color[rgb]{0,0,1}D_{2}}\;\;&{\color[rgb]{0,0.6,0}0}\\\
0&{\color[rgb]{0,0,1}O_{2}^{T}}&{\color[rgb]{0,0,1}D_{3}}\\\
\end{bmatrix}\quad P_{r}=-\begin{bmatrix}0&0&0\\\
{\color[rgb]{0,0.6,0}O_{1}^{T}}&\;\;0\;\;&{\color[rgb]{0,0.6,0}O_{2}}\\\
0&0&0\\\ \end{bmatrix}.\\\ \end{split}$ (11)

The stair is invertible if and only if its diagonal blocks, $D$, are
invertible. If so, it’s inverse has the same stair shaped sparsity pattern
drawing from the same or neighboring block-rows of $S$, enabling efficient
parallel computation [29, 54]:

$\begin{split}\Psi^{-1}&=D^{-1}(2D-\Psi)D^{-1}\\\
\Psi_{l}^{-1}&=\begin{bmatrix}{\color[rgb]{0,0,1}D_{1}^{-1}}&{\color[rgb]{0,0.6,0}0}&0\\\
{\color[rgb]{0,0,1}-D_{2}^{-1}O_{1}^{T}D_{1}^{-1}}&\;\;{\color[rgb]{0,0,1}D_{2}^{-1}}\;\;&{\color[rgb]{0,0,1}-D_{2}^{-1}O_{2}D_{3}^{-1}}\\\
0&{\color[rgb]{0,0.6,0}0}&{\color[rgb]{0,0,1}D_{3}^{-1}}\\\ \end{bmatrix}\\\
\Psi_{r}^{-1}&=\begin{bmatrix}{\color[rgb]{0,0,1}D_{1}^{-1}}&{\color[rgb]{0,0,1}-D_{1}^{-1}O_{1}D_{2}^{-1}}&0\\\
{\color[rgb]{0,0.6,0}0}&\;\;{\color[rgb]{0,0,1}D_{2}^{-1}}\;\;&{\color[rgb]{0,0.6,0}0}\\\
0&{\color[rgb]{0,0,1}-D_{3}^{-1}O_{2}^{T}D_{2}^{-1}}&{\color[rgb]{0,0,1}D_{3}^{-1}}\\\
\end{bmatrix}.\end{split}$ (12)

While, the left and right stair preconditioners are not symmetric, and thus
cannot be used in the context of PCG, the additive polynomial preconditioner
can instead be used:

$\Phi_{add}^{-1}=\frac{1}{2}(\Psi_{l}^{-1}+\Psi_{r}^{-1}).$ (13)

However, as we shall see later on, this will result in a declustered spectrum
that could negatively impact the performance of PCG. In the remainder of this
section we prove a few lemmas and theorems about the eigenpairs of
$\Psi_{l}^{-1}$, $\Psi_{r}^{-1}$, and $\Phi_{add}^{-1}S$ which we will need
for later proofs.222We note that throughout this paper we analyze the spectrum
of $\Psi^{-1}S$ even though the resulting matrix is not guaranteed to be
symmetric positive definite. This is possible because through the Cholesky
factorization, $\Psi=LL^{T}$, we can instead from the symmetric positive
definite $L^{-1}SL^{-T}$, and the matrices $L^{-1}SL^{-T}$ and $\Psi^{-1}S$
are similar.

###  III-A Eigenpairs of $\Psi_{l}$ and $\Psi_{l}$

We first prove a few lemmas about the eigenpairs of the left and right stair
preconditioners.

###### Lemma III.1

Given the stair-splittings of a symmetric block tridiagonal matrix
$S=\Psi_{l}-P_{l}=\Psi_{r}-P_{r}$ for a $n\times n$ block $S$, where each
block is $m\times m$, $\Psi_{l}^{-1}P_{l}$ and $\Psi_{r}^{-1}P_{r}$ have the
same spectrum.

###### Proof:

If $(\lambda,v)$ is an eigenpair of $P_{r}\Psi_{r}^{-1}$ and $\lambda\neq 0$:

$\Psi_{r}^{-1}P_{r}(\Psi_{r}^{-1}v)=\Psi_{r}^{-1}(\lambda v),$ (14)

and $(\lambda,\Psi_{r}^{-1}v)$ is an eigenpair of $\Psi_{r}^{-1}P_{r}$. Since
$P_{r}$ is not invertible, $0$ is an eigenvalue for both $P_{r}\Psi_{r}^{-1}$
and $\Psi_{r}^{-1}P_{r}$. Also as $P_{r}$ and $\Psi_{r}^{-1}$ are equal
dimensioned square matrices then by Theorem 1.3.22 of [55],
$P_{r}\Psi_{r}^{-1}$ and $\Psi_{r}^{-1}P_{r}$ have the same spectrum. Finally,
$(\Psi_{l}^{-1}P_{l})^{T}=P_{r}\Psi_{r}^{-1}$ and every matrix has the same
spectrum as its transpose. ∎

###### Lemma III.2

For
$v^{T}=(v_{1}^{T},v_{2}^{T},\dots,v_{i}^{T},\dots,v_{n}^{T})\in\mathbb{R}^{nm}$
and $v_{i}\in\mathbb{R}^{m}$ for $i=1,2,\dots,n$, we denote
$v_{e}^{T}=(0,v_{2}^{T},\dots,0,v_{2j}^{T},\dots)$ and
$v_{o}^{T}=(v_{1}^{T},0,\dots,v_{2j+1}^{T},0,\dots)$ such that
$v=v_{e}+v_{o}$.

If $(\lambda,v=v_{e}+v_{o})$ is an eigenpair of $\Psi_{l}^{-1}P_{l}$ and
$\lambda\neq 0$ , then $(\lambda,u=v_{e}+\lambda v_{o})$ is an eigenpair of
$\Psi_{r}^{-1}P_{r}$.

If $(\lambda,u=u_{e}+u_{o})$ is an eigenpair of $\Psi_{r}^{-1}P_{r}$ and
$\lambda\neq 0$, then $(\lambda,v=\lambda u_{e}+u_{o})$ is an eigenpair of
$\Psi_{l}^{-1}P_{l}$.

###### Proof:

$\Psi_{l}^{-1}P_{l}v_{o}=0$ for all $v_{o}$ because of the zero columns of
$\Psi_{l}^{-1}P_{l}$ (see Equation 45 in the appendix). Similarly,
$P_{r}v_{e}=0$ and $\Psi_{r}^{-1}P_{r}v_{e}=0$ (see Equation 46 in the
appendix).

Suppose $(\lambda,v=v_{e}+v_{o})$ is an eigenpair of $\Psi_{l}^{-1}P_{l}$,
then

$\Psi_{l}^{-1}P_{l}v_{e}=\Psi_{l}^{-1}P_{l}v=\lambda(v_{e}+v_{o}).$ (15)

If $\lambda\neq 0$, we split Equation 15 into two parts,

$(\Psi_{l}^{-1}-D^{-1})P_{l}v=(\Psi_{l}^{-1}-D^{-1})P_{l}v_{e}=\lambda v_{e}$
(16) $D^{-1}P_{l}v=D^{-1}P_{l}v_{e}=\lambda v_{o}\quad\;\;$ (17)

where $D$ is the diagonal blocks of $S$ (see Equations 47 and 48 in the
appendix for more details).

Since $P_{r}=\Psi_{r}-S=D-\Psi_{l}$, then:

$\begin{split}D^{-1}P_{r}&=I-D^{-1}\Psi_{l}=(\Psi_{l}^{-1}-D^{-1})\Psi_{l}.\end{split}$
(18)

This means that:

$\begin{split}D^{-1}P_{r}v&=(\Psi_{l}^{-1}-D^{-1})\Psi_{l}v\\\
&=(\Psi_{l}^{-1}-D^{-1})\Psi_{l}(\frac{1}{\lambda}\Psi_{l}^{-1}P_{l}v)\\\
&=\frac{1}{\lambda}(\Psi_{l}^{-1}-D^{-1})P_{l}v\\\ &=v_{e}\end{split}$ (19)

Using the inverse of the stair matrix (Equation 12):

$\Psi_{r}^{-1}-D^{-1}=D^{-1}(D-\Psi_{r})D^{-1}=D^{-1}P_{l}D^{-1}$ (20)
$\begin{split}(\Psi_{r}^{-1}-D^{-1})P_{r}v&=D^{-1}P_{l}D^{-1}P_{r}v\\\
&=D^{-1}P_{l}v_{e}\\\ &=\lambda v_{o}\end{split}$ (21)

Combining Equations 19 and 21 as well as make use of the fact $P_{r}v_{e}=0$,
we obtain the following result:

$\begin{split}\Psi_{r}^{-1}P_{r}(v_{e}+\lambda
v_{o})&=\Psi_{r}^{-1}P_{r}(\lambda v)\\\
&=\lambda(D^{-1}P_{r}v+(\Psi_{r}^{-1}-D^{-1})P_{r}v)\\\
&=\lambda(v_{e}+\lambda v_{o})\end{split}$ (22)

The second statement follows from Lemma 4.1 or can be proved with the same
strategy as above. ∎

Finally, we note that the eigenvectors with eigenvalues at $0$ can be divided
into two types. The first type of these are the eigenvectors $v_{o}$ and
$u_{e}$ which are formed because of the zero columns of $\Psi_{l}^{-1}P_{l}$
and $\Psi_{r}^{-1}P_{r}$ (see Equations 45 and 46 in the appendix). The second
type of these are the eigenvectors $v_{e}$ and $u_{o}$ such that
$P_{l}v_{e}=0$ and $P_{r}u_{o}=0$.

If $n$ is even, then $\Psi_{l}^{-1}P_{l}$ and $\Psi_{r}^{-1}P_{r}$ have
exactly $m\frac{n}{2}$ eigenpairs of the first type. If $n$ is odd,
$\Psi_{l}^{-1}P_{l}$ will have $m\frac{n+1}{2}$ while $\Psi_{r}^{-1}P_{r}$
will have $m$ less such eigenpair and, instead, have $m$ more eigenpairs in
the form of $(0,u_{o})$ since, in this case, $P_{r}u_{o}=0$, which has non-
zero solutions even when all $O$s are invertible:

$\left\\{\begin{array}[]{lr}O_{1}^{T}u_{1}+O_{2}u_{3}=0\\\
O_{3}^{T}u_{3}+O_{4}u_{5}=0\\\ \cdots\\\
O_{n-2}^{T}u_{n-2}+O_{n-1}u_{n}=0\end{array}\right.$ (23)

In the case of most integrators, as discussed before, $O$s are invertible, and
$P_{l}v_{e}=0$ if and only if $v_{e}=0$. Hence, there are exactly
$m\lceil\frac{n}{2}\rceil$ eigenpairs with eigenvalues at $0$ and
$m\lfloor\frac{n}{2}\rfloor$ eigenpairs with non-zero eigenvalues.

###  III-B Eigenpairs of $\Phi_{add}$

We next construct the $mn$ eigenpairs of $\Phi_{add}^{-1}S$ from those of
$\Psi_{l}^{-1}P_{l}$ and $\Psi_{r}^{-1}P_{r}$.

###### Theorem III.3

If $(0,v_{e})$ is an eigenpair of $\Psi_{l}^{-1}P_{l}$ then $(1,v_{e})$ is an
eigenpair of $\Phi_{add}^{-1}S$ and if $(0,u_{o})$ is an eigenpair of
$\Psi_{r}^{-1}P_{r}$ then $(1,u_{o})$ are eigenpairs of $\Phi_{add}^{-1}S$.

###### Proof:

Let $(0,v_{e})$ be an eigenpair of $\Psi_{l}^{-1}P_{l}$, then:

$\Psi_{l}^{-1}Sv_{e}=\Psi_{l}^{-1}(\Psi_{l}-P_{l})v_{e}=v_{e}-0v_{e}=v_{e}$
(24) $\Psi_{r}^{-1}Sv_{e}=\Psi_{r}^{-1}(\Psi_{r}-P_{r})v_{e}=v_{e}-0=v_{e}.$
(25)
$\Phi_{add}^{-1}Sv_{e}=\frac{1}{2}(\Psi_{l}^{-1}+\Psi_{r}^{-1})Sv_{e}=\frac{1}{2}(v_{e}+v_{e})=v_{e}.$
(26)

Let $(0,u_{o})$ be an eigenpair of $\Psi_{r}^{-1}P_{r}$, then

$\Psi_{l}^{-1}Su_{o}=\Psi_{l}^{-1}(\Psi_{l}-P_{l})u_{o}=u_{o}-0=u_{o}$ (27)
$\Psi_{r}^{-1}Su_{o}=\Psi_{r}^{-1}(\Psi_{r}-P_{r})u_{o}=u_{o}-0u_{o}=u_{o}.$
(28) $\Phi_{add}^{-1}Su_{o}=u_{o}.$ (29)

∎

###### Theorem III.4

Let $(\lambda,v=v_{e}+v_{o})$ be an eigenpair of $\Psi_{l}^{-1}P_{l}$. If
$\lambda\neq 0$, then
$(1-\frac{1}{2}(\lambda\pm\sqrt{\lambda}),v_{e}\pm\sqrt{\lambda}v_{o})$ are
eigenpairs of $\Phi_{add}^{-1}S$.

###### Proof:

Let $(\lambda,v=v_{e}+v_{o})$ be an eigenpair of $\Psi_{l}^{-1}P_{l}$: Since
$\Psi_{l}^{-1}P_{l}v_{o}=0$,

$\Psi_{l}^{-1}P_{l}v_{e}=\lambda(v_{e}+v_{o}).$ (30)

Similarly, Equation 22 and $\Psi_{r}^{-1}P_{r}v_{e}=0$ suggests:

$\Psi_{r}^{-1}P_{r}v_{o}=v_{e}+\lambda v_{o}.$ (31)

We can therefore make an educated guess that the eigenvalues of
$\Phi_{add}^{-1}S=\frac{1}{2}(\Psi_{l}^{-1}+\Psi_{r}^{-1})S=\frac{1}{2}(\Psi_{l}^{-1}S+\Psi_{r}^{-1}S)$
should be of the form $u=av_{e}+bv_{o}$. Following from Equations 30 and 31,
we can therefore say that:

$\begin{split}\Phi_{add}^{-1}Su&=\frac{1}{2}(\Psi_{l}^{-1}S(av_{e}+bv_{o})+\Psi_{r}^{-1}S(av_{e}+bv_{o}))\\\
&=u-\frac{1}{2}\left[b\Psi_{l}^{-1}P_{l}v_{e}+a\Psi_{r}^{-1}P_{r}v_{o}\right]\\\
&=u-\frac{1}{2}\left[b\lambda(v_{e}+v_{o})+a(v_{e}+\lambda v_{o})\right]\\\
&=u-\frac{1}{2}\left[(a+\lambda b)v_{e}+\lambda(a+b)v_{o}\right]\\\
\end{split}$

(32)

Given our assumed form $u=av_{e}+bv_{o}$, the above equation should simplify
to: $\beta av_{e}+\beta bv_{o}$. Thus, the eigenpair $(\lambda,u)$ must
satisfy the linear system:

$\frac{1}{2}\left(\begin{matrix}1&\lambda\\\
\lambda&\lambda\end{matrix}\right)\left(\begin{matrix}a\\\
b\end{matrix}\right)=(1-\beta)\left(\begin{matrix}a\\\ b\end{matrix}\right)$
(33)

This means we can conclude that:

$\beta=1-\frac{1}{2}(\lambda\pm\sqrt{\lambda}),\quad\quad
b=\pm\sqrt{\lambda}a$ (34)

∎

###### Theorem III.5

If $S$ is symmetric block tridiagonal and positive definite, then
$\Phi_{add}^{-1}S$ has real positive eigenvalues.

###### Proof:

If $S$ is positive definite, then its diagonal block matrix, $D$, is positive
definite and thus:

$\begin{split}v^{T}Sv&>0\\\ v^{T}(\Psi_{l}+\Psi_{r}-D)v&>0\\\
v^{T}(\Psi_{l}+\Psi_{r})v&>0\end{split}$ (35)

Since $\Psi_{l}=\Psi_{r}^{T}$ we can also note that $\Psi_{l}$, $\Psi_{r}$,
$\Psi_{l}^{-1}$, and $\Psi_{r}^{-1}$ are all positive definite as:

$v^{T}(\Psi_{l}+\Psi_{r})v=2v^{T}\Psi_{l}v=2v^{T}\Psi_{r}v$ (36)

As $(\Psi_{l}^{-1})^{T}=\Psi_{r}^{-1}$,
$\Phi_{add}^{-1}=\frac{1}{2}(\Psi_{l}^{-1}+\Psi_{r}^{-1})$ is symmetric and
positive definite. $\Phi_{add}^{-1}S$ is then similar to the symmetric
positive definite matrix
$\frac{1}{2}(\Psi_{l}^{-1}+\Psi_{r}^{-1})^{1/2}S(\Psi_{l}^{-1}+\Psi_{r}^{-1})^{1/2}$
and must also have positive real eigenvalues. We note that this is a special
case of Theorem 3.4 of [34] and Theorem 2.2 and Corollary 2.3 of [46]. ∎

##  IV The Symmetric Stair Precontioner

In this section, we introduce a new symmetric parallel preconditioner for
symmetric positive definite block tridiagonal matrices which improves upon
existing stair based preconditoners, the _symmetric stair preconditioner_ :

$\Phi^{-1}_{sym}=\Psi_{l}^{-1}+\Psi_{r}^{-1}-D^{-1}.$ (37)

This new preconditioner can also simply be thought of as taking
$\Psi_{l}^{-1}$ or $\Psi_{r}^{-1}$ and copying the off diagonal blocks across
the diagonal and takes the form:

$\Phi^{-1}_{sym}=\begin{bmatrix}{\color[rgb]{0,0,1}D_{1}^{-1}}&{\color[rgb]{0,0,1}-D_{1}^{-1}O_{1}D_{2}^{-1}}&0\\\
{\color[rgb]{0,0,1}-D_{2}^{-1}O_{1}^{T}D_{1}^{-1}}&\;\;{\color[rgb]{0,0,1}D_{2}^{-1}}\;\;&{\color[rgb]{0,0,1}-D_{2}^{-1}O_{2}D_{3}^{-1}}\\\
0&{\color[rgb]{0,0,1}-D_{3}^{-1}O_{2}^{T}D_{2}^{-1}}&{\color[rgb]{0,0,1}D_{3}^{-1}}\\\
\end{bmatrix}\\\ $

(38)

In the remainder of this section, we will prove that $\Phi_{sym}^{-1}S$ is
admissible for use with PCG, and that $\Phi_{sym}^{-1}S$ has a more clustered
spectrum than $\Phi_{add}^{-1}S$.

We begin by showing that $\Phi_{sym}^{-1}S$ and $\Psi_{l}^{-1}S$ (or
$\Psi_{r}^{-1}S$) share the same eigenvalues. And, that the multiplicities of
the eigenvalues of $\Phi_{sym}^{-1}S$ that are less than one are doubled as
compared to $\Psi_{l}^{-1}S$ (or $\Psi_{r}^{-1}S$).

###### Theorem IV.1

If $(0,v_{e})$ is an eigenpair of $\Psi_{l}^{-1}P_{l}$ then $(1,v_{e})$ is an
eigenpair of $\Phi_{sym}^{-1}S$. If $(0,u_{o})$ is an eigenpair of
$\Psi_{r}^{-1}P_{r}$ then $(1,u_{o})$ are eigenpairs of $\Phi_{sym}^{-1}S$.

###### Proof:

Let $(0,v_{e})$ be an eigenpair of $\Psi_{l}^{-1}P_{l}$, then by Equations 25
and 48 in the appendix:

$\begin{split}\Phi_{sym}^{-1}Sv_{e}&=(\Psi_{r}^{-1}+(\Psi_{l}^{-1}-D^{-1}))Sv_{e}\\\
&=v_{e}+(\Psi_{l}^{-1}-D^{-1})(\Psi_{l}-P_{l})v_{e}\\\
&=v_{e}-(\Psi_{l}^{-1}-D^{-1})P_{l}v_{e}+(\Psi_{l}^{-1}-D^{-1})\Psi_{l}v_{e}\\\
&=v_{e}-0v_{e}+0=v_{e}\end{split}$

(39)

Similarly, let $(0,u_{o})$ be an eigenpair of $\Psi_{r}^{-1}P_{r}$, then:

$\begin{split}\Phi_{sym}^{-1}Su_{o}&=(\Psi_{l}^{-1}+(\Psi_{r}^{-1}-D^{-1}))Su_{o}\\\
&=u_{o}+(\Psi_{r}^{-1}-D^{-1})(\Psi_{r}-P_{r})u_{o}\\\
&=u_{o}-(\Psi_{r}^{-1}-D^{-1})P_{r}u_{o}+(\Psi_{l}^{-1}-D^{-1})\Psi_{r}u_{o}\\\
&=u_{o}-0u_{o}+0=u_{o}\end{split}$

(40)

∎

###### Theorem IV.2

If $(\lambda,v=v_{e}+v_{o})$ is an eigenpair of $\Psi_{l}^{-1}P_{l}$, then
$(1-\lambda,v_{o})$ and $(1-\lambda,v_{e})$ are eigenpairs of
$\Phi_{sym}^{-1}S$.

###### Proof:

Using the splittings in Equations 16 and 17:

$\begin{split}\Phi_{sym}^{-1}Sv_{e}&=(\Psi_{r}^{-1}+(\Psi_{l}^{-1}-D^{-1}))Sv_{e}\\\
&=\Psi_{r}^{-1}(\Psi_{r}-P_{r})v_{e}+(\Psi_{l}^{-1}-D^{-1})(\Psi_{l}-P_{l})v_{e}\\\
&=v_{e}+(\Psi_{l}^{-1}-D^{-1})\Psi_{l}v_{e}-(\Psi_{l}^{-1}-D^{-1})P_{l}v_{e}\\\
&=v_{e}+D^{-1}P_{r}v_{e}-\lambda v_{e}\\\ &=(1-\lambda)v_{e}.\end{split}$

(41)

Hence, $(1-\lambda,v_{e})$ is an eigenpair of $\Phi_{sym}^{-1}S$. The proof
for $(1-\lambda,v_{o})$ follows the same structure. ∎

Finally, we will prove that $\Phi_{sym}^{-1}S$ is admissible for PCG, and that
$\Phi_{sym}^{-1}S$ has a smaller condition number than $\Phi_{add}^{-1}S$.

###### Theorem IV.3

If $S$ is symmetric block tridiagonal and positive definite, then:

  1. 1.

$\Phi_{sym}^{-1}S$ has real eigenvalues.

  2. 2.

$\Phi_{sym}^{-1}P$ has non-negative real eigenvalues and
$\rho(\Phi_{sym}^{-1}P)<1$.

  3. 3.

$\Phi_{sym}^{-1}S$ is positive definite and
$\lambda(\Phi_{sym}^{-1}S)\in(0,1]$.

$\lambda_{\max}(\Phi_{sym}^{-1}S)=1-\lambda_{\min}(\Phi_{sym}^{-1}P)\leq
1<1-\frac{1}{2}\left(\lambda^{+}(\Psi^{-1}P)-\sqrt{\lambda^{+}(\Psi^{-1}P)}\right)\leq\lambda_{\max}(\Phi_{add}^{-1}S)\leq\frac{9}{8}$
(42)

$\lambda_{\min}(\Phi_{sym}^{-1}S)=1-\lambda^{+}_{\max}(\Phi_{sym}^{-1}P)>1-\frac{1}{2}\left(\lambda^{+}_{\max}(\Psi^{-1}P)+\sqrt{\lambda^{+}_{\max}(\Psi^{-1}P)}\right)=\lambda_{\min}(\Phi_{add}^{-1}S)>0$
(43)

###### Proof:

We prove each part in order:

  1. 1.

As $S$ is positive definite, then $\Phi_{sym}^{-1}S$ is similar to the
symmetric matrix $S^{1/2}\Phi_{sym}^{-1}S^{1/2}$. Therefore, the eigenvalues
of $\Phi_{sym}^{-1}S$ are real.

  2. 2.

As $\lambda\in\mathbb{R}$ from point one, then as we again derive $\lambda$
from $\Psi$, by Equation 34, $1-\frac{1}{2}(\lambda\pm\sqrt{\lambda})$ will be
real and positive if and only if $0\leq\lambda<1$. Thus, $\Phi_{sym}^{-1}P$
has non-negative real eigenvalues and $\rho(\Phi_{sym}^{-1}P)<1$. This also
shows that the stair splitting is convergent for a symmetric positive definite
$S$.

  3. 3.

From points one and two and Theorem 4.1 and Theorem 4.2, $1-\lambda\in(0,1]$
and $\Phi_{sym}^{-1}S$ is positive definite.

∎

As a result of Theorem 3.4, Theorem 3.5, and Theorem 4.3, we have the
following relationships between the maximum (Equation 42) and minimum
(Equation 43) eigenvalues of the two preconditioners, where
$\lambda^{+}_{\max}(\Psi_{x}^{-1}P)$ and $\lambda^{+}_{\min}(\Psi_{x}^{-1}P)$
denotes the largest and smallest non-zero eigenvalues of $\Psi_{x}$. As we can
see from Equations 42 and 43, the symmetric stair preconditioner results in a
smaller condition number as the eigenvalues are of the range:

$\begin{split}\lambda(\Phi_{sym}^{-1}S)&\in(0,1],\\\
\lambda(\Phi_{add}^{-1}S)&\in\left(0,\frac{9}{8}\right].\end{split}$ (44)

##  V Numerical Results

In this section we present a numerical evaluation of our symmetric stair
preconditioner and compare it to the additive stair preconditioner and other
parallel preconditioners from the literature on representative trajectory
optimization tasks. We evaluate the spread of the eigenvalues both in absolute
terms and with regard to the the relative condition number after
preconditioning, as well as the resulting number of iterations of PCG needed
for convergence. We use the canonical pendulum and cart pole swing up problems
as well as a problem to compute a motion across its workspace for a Kuka LBR
IIWA-14 manipulator. Source code accompanying this evaluation, including all
hyperparameter values used in these experiments, can be found at
github.com/a2r-lab/SymStair.

Figure 1: Distribution of the Eigenvalues of the additive and symmetric stair
preconditioners matching the theoretical results in Equation 44.

Figure 1 shows the distribution of the eigenvalues of the additive and
symmetric stair preconditioners along with the line for 0.0 highlighted in
black, 1.0 in green, and $\frac{9}{8}=1.125$ in purple. We see that across all
problems our numerical results match the theoretical results in Equation 44.

This theoretical result directly leads to a better relative condition number
for our preconditioner as compared to alternative parallel preconditioners
from the literature. As shown in Figure 2, when normalizing the condition
number to the Jacobi preconditioner equaling 1 for each system, to enable
comparison across the three problems which have large differences in the
absolute condition number, we see that the symmetric stair preconditioner is
the most performant. In fact, it is not only able to reduce the relative
condition number by 76-89% as compared to the standard Jacobi preconditioner,
but also outperforms the best parallel preconditioner from the literature, the
additive stair preconditioner, by 33-34%.

This improvement in numerical conditioning also drastically reduces the number
of PCG iterations needed for the problem to converge, enabling faster linear
system solves. This is crucial for real-time trajectory optimization as each
trajectory optimization solve requires the solution of many linear systems. In
particular, as shown in Figure 3, when solving the first KKT system for our
three target trajectory optimization problems, PCG using the symmetric stair
preconditioner requires 51-68% fewer iterations than the Jacobi
preconditioner, and 17-25% fewer iterations than the next best parallel
preconditioner to converge to a solution under the same exit tolerance.

$\Psi_{l}^{-1}P_{l}=\displaystyle\left(\begin{matrix}0&-D_{1}^{-1}O_{1}&0&0&0&0\\\
0&D_{2}^{-1}O_{2}D_{3}^{-1}O_{2}^{T}+D_{2}^{-1}O_{1}^{T}D_{1}^{-1}O_{1}&0&D_{2}^{-1}O_{2}D_{3}^{-1}O_{3}&0&0\\\
0&-D_{3}^{-1}O_{2}^{T}&0&-D_{3}^{-1}O_{3}&0&0\\\
0&D_{4}^{-1}O_{3}^{T}D_{3}^{-1}O_{2}^{T}&0&D_{4}^{-1}O_{4}D_{5}^{-1}O_{4}^{T}+D_{4}^{-1}O_{3}^{T}D_{3}^{-1}O_{3}&0&D_{4}^{-1}O_{4}D_{5}^{-1}O_{5}\\\
0&0&0&-D_{5}^{-1}O_{4}^{T}&0&-D_{5}^{-1}O_{5}\\\
0&0&0&D_{6}^{-1}O_{5}^{T}D_{5}^{-1}O_{4}^{T}&0&D_{6}^{-1}O_{5}^{T}D_{5}^{-1}O_{5}\end{matrix}\right)$

(45)

$\Psi_{r}^{-1}P_{r}=\displaystyle\left(\begin{matrix}D_{1}^{-1}O_{1}D_{2}^{-1}O_{1}^{T}&0&D_{1}^{-1}O_{1}D_{2}^{-1}O_{2}&0&0&0\\\
-D_{2}^{-1}O_{1}^{T}&0&-D_{2}^{-1}O_{2}&0&0&0\\\
D_{3}^{-1}O_{2}^{T}D_{2}^{-1}O_{1}^{T}&0&D_{3}^{-1}O_{3}D_{4}^{-1}O_{3}^{T}+D_{3}^{-1}O_{2}^{T}D_{2}^{-1}O_{2}&0&D_{3}^{-1}O_{3}D_{4}^{-1}O_{4}&0\\\
0&0&-D_{4}^{-1}O_{3}^{T}&0&-D_{4}^{-1}O_{4}&0\\\
0&0&D_{5}^{-1}O_{4}^{T}D_{4}^{-1}O_{3}^{T}&0&D_{5}^{-1}O_{5}D_{6}^{-1}O_{5}^{T}+D_{5}^{-1}O_{4}^{T}D_{4}^{-1}O_{4}&0\\\
0&0&0&0&-D_{6}^{-1}O_{5}^{T}&0\end{matrix}\right)$

(46)

$(\Psi_{l}^{-1}-D^{-1})P_{l}=\displaystyle\left(\begin{matrix}0&0&0&0&0&0\\\
0&D_{2}^{-1}O_{2}D_{3}^{-1}O_{2}^{T}+D_{2}^{-1}O_{1}^{T}D_{1}^{-1}O_{1}&0&D_{2}^{-1}O_{2}D_{3}^{-1}O_{3}&0&0\\\
0&0&0&0&0&0\\\
0&D_{4}^{-1}O_{3}^{T}D_{3}^{-1}O_{2}^{T}&0&D_{4}^{-1}O_{4}D_{5}^{-1}O_{4}^{T}+D_{4}^{-1}O_{3}^{T}D_{3}^{-1}O_{3}&0&D_{4}^{-1}O_{4}D_{5}^{-1}O_{5}\\\
0&0&0&0&0&0\\\
0&0&0&D_{6}^{-1}O_{5}^{T}D_{5}^{-1}O_{4}^{T}&0&D_{6}^{-1}O_{5}^{T}D_{5}^{-1}O_{5}\end{matrix}\right)$

(47)

$D^{-1}P_{l}=\displaystyle\left(\begin{matrix}0&-D_{1}^{-1}O_{1}&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&-D_{3}^{-1}O_{2}^{T}&0&-D_{3}^{-1}O_{3}&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&-D_{5}^{-1}O_{4}^{T}&0&-D_{5}^{-1}O_{5}\\\
0&0&0&0&0&0\end{matrix}\right)\quad\quad(\Psi_{l}^{-1}-D^{-1})\Psi_{l}=\displaystyle\left(\begin{matrix}0&0&0&0&0&0\\\
-D_{2}^{-1}O_{1}^{T}&0&-D_{2}^{-1}O_{2}&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&-D_{4}^{-1}O_{3}^{T}&0&-D_{4}^{-1}O_{4}&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&-D_{6}^{-1}O_{5}^{T}&0\end{matrix}\right)$

(48)

Figure 2: The relative condition number resulting from different
preconditioners, normalized to the value of the Jacobi preconditioner, showing
the improved performance of the symmetric stair preconditioner. Figure 3: The
number of PCG iterations required for convergence to the same exit tolerance
across different preconditioners and problems, again showing the improved
performance of the symmetric stair preconditioner.

##  VI Conclusion and Future Work

In this work we present a new parallel-friendly symmetric stair
preconditioner. Through both proofs and numerical experiments, we show that
our preconditioner has advantageous theoretical and practical properties that
enable fast iterative linear system solves. Across three benchmark tasks, our
preconditoner provides up to a 34% reduction in condition number and 25%
reduction in PCG iterations.

In future work we hope to explore the theoretical and practical performance
tradeoffs of higher order polynomial preconditioners on heterogeneous
hardware, as well as extend our analysis to provide bounds on the number of
PCG iterations under various exit tolerances and preconditioners.

##  VII Appendix

Some useful results revealing the sparsity patterns for the $6\times 6$ block
cases are listed in Equations 45\- 48.

## References

  * [1] J. T. Betts, _Practical Methods for Optimal Control Using Nonlinear Programming_ , ser. Advances in Design and Control. Society for Industrial and Applied Mathematics (SIAM), vol. 3.
  * [2] Y. Tassa, T. Erez, and E. Todorov, “Synthesis and Stabilization of Complex Behaviors through Online Trajectory Optimization,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_.
  * [3] W. Xi and C. D. Remy, “Optimal Gaits and Motions for Legged Robots,” in _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_.
  * [4] M. Posa, S. Kuindersma, and R. Tedrake, “Optimization and stabilization of trajectories for constrained dynamical systems,” in _Proceedings of the International Conference on Robotics and Automation (ICRA)_. IEEE, pp. 1366–1373.
  * [5] T. Apgar, P. Clary, K. Green, A. Fern, and J. Hurst, “Fast Online Trajectory Optimization for the Bipedal Robot Cassie,” in _Proceedings of Robotics: Science and Systems_.
  * [6] T. Howell, B. Jackson, and Z. Manchester, “Altro: A fast solver for constrained trajectory optimization,” in _Proceedings of (IROS) IEEE/RSJ International Conference on Intelligent Robots and Systems_ , November 2019, pp. 7674 – 7679.
  * [7] S. Kuindersma, “Taskable agility: Making useful dynamic behavior easier to create,” Princeton Robotics Seminar, April 2023.
  * [8] J. T. Betts and W. P. Huffman, “Trajectory optimization on a parallel processor,” vol. 14, no. 2, pp. 431–439.
  * [9] D. Kouzoupis, R. Quirynen, B. Houska, and M. Diehl, “A Block Based ALADIN Scheme for Highly Parallelizable Direct Optimal Control,” in _Proceedings of the American Control Conference_.
  * [10] M. Giftthaler, M. Neunert, M. Stäuble, J. Buchli, and M. Diehl, “A family of iterative gauss-newton shooting methods for nonlinear optimal control,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 1–9.
  * [11] F. Farshidian, E. Jelavic, A. Satapathy, M. Giftthaler, and J. Buchli, “Real-time motion planning of legged robots: A model predictive control approach,” in _2017 IEEE-RAS 17th International Conference on Humanoid Robotics_.
  * [12] A. Astudillo, J. Gillis, G. Pipeleers, W. Decré, and J. Swevers, “Speed-up of nonlinear model predictive control for robot manipulators using task and data parallelism,” in _2022 IEEE 17th International Conference on Advanced Motion Control (AMC)_ , 2022, pp. 201–206.
  * [13] T. Antony and M. J. Grant, “Rapid Indirect Trajectory Optimization on Highly Parallel Computing Architectures,” vol. 54, no. 5, pp. 1081–1091.
  * [14] B. Plancher and S. Kuindersma, “A Performance Analysis of Parallel Differential Dynamic Programming on a GPU,” in _International Workshop on the Algorithmic Foundations of Robotics (WAFR)_.
  * [15] ——, “Realtime model predictive control using parallel ddp on a gpu,” in _Toward Online Optimal Control of Dynamic Robots Workshop at the 2019 International Conference on Robotics and Automation (ICRA)_ , Montreal, Canada, May. 2019.
  * [16] Z. Pan, B. Ren, and D. Manocha, “Gpu-based contact-aware trajectory optimization using a smooth force model,” in _Proceedings of the 18th Annual ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , ser. SCA ’19. New York, NY, USA: ACM, 2019, pp. 4:1–4:12.
  * [17] B. Plancher, S. M. Neuman, T. Bourgeat, S. Kuindersma, S. Devadas, and V. J. Reddi, “Accelerating robot dynamics gradients on a cpu, gpu, and fpga,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2335–2342, 2021\.
  * [18] B. Plancher, S. M. Neuman, R. Ghosal, S. Kuindersma, and V. J. Reddi, “GRiD: GPU-Accelerated Rigid Body Dynamics with Analytical Gradients,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, pp. 6253–6260. [Online]. Available: https://ieeexplore.ieee.org/document/9812384/
  * [19] J. Sacks, D. Mahajan, R. C. Lawson, and H. Esmaeilzadeh, “Robox: an end-to-end solution to accelerate autonomous control in robotics,” in _Proceedings of the 45th Annual International Symposium on Computer Architecture_. IEEE Press, 2018, pp. 479–490.
  * [20] S. M. Neuman, B. Plancher, T. Bourgeat, T. Tambe, S. Devadas, and V. J. Reddi, “Robomorphic computing: A design methodology for domain-specific accelerators parameterized by robot morphology,” ser. ASPLOS 2021. New York, NY, USA: Association for Computing Machinery, 2021, p. 674–686. [Online]. Available: https://doi.org/10.1145/3445814.3446746
  * [21] S. M. Neuman, R. Ghosal, T. Bourgeat, B. Plancher, and V. J. Reddi, “Roboshape: Using topology patterns to scalably and flexibly deploy accelerators across robots,” in _Proceedings of the 50th Annual International Symposium on Computer Architecture_ , ser. ISCA ’23. New York, NY, USA: Association for Computing Machinery, 2023. [Online]. Available: https://doi.org/10.1145/3579371.3589104
  * [22] Y. Yang, X. Chen, and Y. Han, “Rbdcore: Robot rigid body dynamics accelerator with multifunctional pipelines,” _arXiv preprint arXiv:2307.02274_ , 2023\.
  * [23] H. Esmaeilzadeh, E. Blem, R. St. Amant, K. Sankaralingam, and D. Burger, “Dark Silicon and the End of Multicore Scaling,” in _Proceedings of the 38th Annual International Symposium on Computer Architecture_ , ser. ISCA ’11. ACM, pp. 365–376.
  * [24] G. Venkatesh, J. Sampson, N. Goulding, S. Garcia, V. Bryksin, J. Lugo-Martinez, S. Swanson, and M. B. Taylor, “Conservation Cores: Reducing the Energy of Mature Computations,” in _Proceedings of the Fifteenth Edition of ASPLOS on Architectural Support for Programming Languages and Operating Systems_ , ser. ASPLOS XV. ACM, pp. 205–218.
  * [25] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “Osqp: An operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020.
  * [26] G. Frison and M. Diehl, “Hpipm: a high-performance quadratic programming framework for model predictive control,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 6563–6569, 2020.
  * [27] S. C. Eisenstat, “Efficient implementation of a class of preconditioned conjugate gradient methods,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 2, no. 1, pp. 1–4, 1981.
  * [28] Y. Saad, _Iterative Methods for Sparse Linear Systems_. Society for Industrial and Applied Mathematics, Apr. 2003.
  * [29] B. Plancher, “GPU Acceleration for Real-Time, Whole-Body, Nonlinear Model Predictive Controla,” Ph.D. dissertation, Harvard University, United States – Massachusetts, 2022.
  * [30] M. Schubiger, G. Banjac, and J. Lygeros, “GPU acceleration of ADMM for large-scale quadratic programming,” _Journal of Parallel and Distributed Computing_ , vol. 144, pp. 55–67, Oct. 2020.
  * [31] R. Helfenstein and J. Koko, “Parallel preconditioned conjugate gradient algorithm on gpu,” _Journal of Computational and Applied Mathematics_ , vol. 236, no. 15, pp. 3584–3590, 2012.
  * [32] J. Nocedal and S. J. Wright, _Numerical Optimization_ , 2nd ed. Springer.
  * [33] J. R. Shewchuk, “An introduction to the conjugate gradient method without the agonizing pain,” USA, Tech. Rep., 1994.
  * [34] H.-B. Li, T.-Z. Huang, Y. Zhang, X.-P. Liu, and H. Li, “On some new approximate factorization methods for block tridiagonal matrices suitable for vector and parallel processors,” _Mathematics and Computers in Simulation_ , vol. 79, no. 7, pp. 2135–2147, Mar. 2009.
  * [35] W. Karush, “Minima of functions of several variables with inequalities as side conditions,” Master’s thesis, Department of Mathematics, University of Chicago, Chicago, IL, USA, 1939.
  * [36] H. W. Kuhn and A. W. Tucker, “Nonlinear programming,” in _Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1950_. Berkeley and Los Angeles: University of California Press, 1951, pp. 481–492.
  * [37] H. K. Khalil, _Nonlinear Systems_. Prentice Hall, 2002.
  * [38] J. Bolz, I. Farmer, E. Grinspun, and P. Schröoder, “Sparse Matrix Solvers on the GPU: Conjugate Gradients and Multigrid,” in _ACM SIGGRAPH 2003 Papers_ , ser. SIGGRAPH ’03. ACM, pp. 917–924. [Online]. Available: http://doi.acm.org/10.1145/1201775.882364
  * [39] H. Liu, J.-H. Seo, R. Mittal, and H. H. Huang, “Gpu-accelerated scalable solver for banded linear systems,” in _2013 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2013, pp. 1–8.
  * [40] H. Anzt, M. Gates, J. Dongarra, M. Kreutzer, G. Wellein, and M. Köhler, “Preconditioned krylov solvers on gpus,” _Parallel Computing_ , 05 2017\.
  * [41] H. Anzt, M. Kreutzer, E. Ponce, G. D. Peterson, G. Wellein, and J. Dongarra, “Optimization and performance evaluation of the idr iterative krylov solver on gpus,” _The International Journal of High Performance Computing Applications_ , vol. 32, no. 2, pp. 220–230, 2018.
  * [42] G. Flegar _et al._ , “Sparse linear system solvers on gpus: Parallel preconditioning, workload balancing, and communication reduction,” Ph.D. dissertation, Universitat Jaume I, 2019.
  * [43] M. Schubiger, G. Banjac, and J. Lygeros, “Gpu acceleration of admm for large-scale quadratic programming,” _Journal of Parallel and Distributed Computing_ , vol. 144, pp. 55–67, 2020.
  * [44] J. W. Pearson and J. Pestana, “Preconditioners for krylov subspace methods: An overview,” _GAMM-Mitteilungen_ , vol. 43, no. 4, p. e202000015, 2020.
  * [45] L. Adams, “$m$-Step Preconditioned Conjugate Gradient Methods,” _Society for Industrial and Applied Mathematics. SIAM Journal on Scientific and Statistical Computing_ , vol. 6, no. 2, p. 12, Apr. 1985.
  * [46] L. M. Adams and E. G. Ong, “Additive polynomial preconditioners for parallel computers,” _Parallel Computing_ , vol. 9, no. 3, pp. 333–345, Feb. 1989\.
  * [47] P. Concus, G. H. Golub, and D. P. O’Leary, “Numerical solution of nonlinear elliptic partial differential equations by a generalized conjugate gradient method,” _Computing_ , vol. 19, no. 4, pp. 321–339, Dec. 1978.
  * [48] K. Muzhinji, “Optimal Block Preconditioner for an Efficient Numerical Solution of the Elliptic Optimal Control Problems Using Gmres Solver,” _International Journal of Numerical Analysis and Modeling_ , vol. 20, no. 1, pp. 47–66, June 2023.
  * [49] M. Benzi and A. J. Wathen, “Some Preconditioning Techniques for Saddle Point Problems,” in _Model Order Reduction: Theory, Research Aspects and Applications_ , ser. Mathematics in Industry, W. H. A. Schilders, H. A. van der Vorst, and J. Rommes, Eds. Berlin, Heidelberg: Springer, 2008, pp. 195–211.
  * [50] E. Galligani and V. Ruggiero, “A polynomial preconditioner for block tridiagonal matrices,” _PARALLEL ALGORITHM AND APPLICATIONS_ , vol. 3, no. 3-4, pp. 227–237, 1994.
  * [51] S. Demko, W. F. Moss, and P. W. Smith, “Decay rates for inverses of band matrices,” _Mathematics of computation_ , vol. 43, no. 168, pp. 491–499, 1984.
  * [52] H.-B. Li, T.-Z. Huang, Y. Zhang, X.-P. Liu, and H. Li, “On some new approximate factorization methods for block tridiagonal matrices suitable for vector and parallel processors,” _Mathematics and Computers in Simulation_ , vol. 79, no. 7, pp. 2135–2147, 2009. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378475408003881
  * [53] H.-B. Li, T.-Z. Huang, Y. Zhang, X.-P. Liu, and T.-X. Gu, “Chebyshev-type methods and preconditioning techniques,” _Applied Mathematics and Computation_ , vol. 218, no. 2, pp. 260–270, 2011.
  * [54] H. Lu, “Stair Matrices and Their Generalizations with Applications to Iterative Methods I: A Generalization of the Successive Overrelaxation Method,” vol. 37, no. 1, pp. 1–17. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/S0036142998343294
  * [55] R. A. Horn and C. R. Johnson, _Matrix analysis_. Cambridge university press, 2012.


# Decomposition of the Leinster-Cobbold Diversity Index

Bingzhang Chen, Michael Grinfeld  Department of Mathematics and Statistics,
University of Strathclyde, 26 Richmond Street, Glasgow G1 1XH, UK
bingzhang.chen<EMAIL_ADDRESS>
###### Abstract.

The Leinster and Cobbold diversity index (Leinster and Cobbold, 2012)
possesses a number of merits; in particular, it generalises many existing
indices and defines an effective number. We present a scheme to quantify the
contribution of richness, evenness, and taxonomic similarity to this index.
Compared to the work of (van Dam, 2019), our approach gives unbiased estimates
of both evenness and similarity in a non-homogeneous community. We also
introduce a notion of taxonomic tree equilibration which should be of use in
the description of community structure.

Keywords: diversity, similarity, evenness, richness, decomposition,
Leinster–Cobbold diversity index.

##  1\. Introduction

Measuring biodiversity is a difficult task due to sampling issues and
accounting for missing data, but also as there is no one universally accepted
definition of what biodiversity is (Daly et al., 2018). In ecological
practice, definitions of biodiversity can include contributions from multiple
channels of information such as the number of species (“richness”), dominance
or rarity relations among the constituent species (“evenness”), and measures
of “similarity” among the species (estimated either from taxonomic or
phylogenetic relationships, or from functional traits relationships) (Purvis
and Hector, 2000; Leinster and Cobbold, 2012). Biogeographic patterns of
diversity can depend on the definitions used. For example, (Stuart-Smith et
al., 2013) showed that biodiversity hotspots can shift from the tropics to
higher latitudes if one only considers abundances or also takes account of
functional traits similarity of species.

An outstanding challenge in the field of conservation ecology is to relate the
various aspects of biodiversity data to the functioning of ecosystems
(Maureaud et al., 2019; Hillebrand et al., 2018).

Thus the goal is to construct a biodiversity index that would carry
information about as many aspects of diversity as possible. This goal has been
actively pursued (Rao, 1982; Leinster and Cobbold, 2012; Chao et al., 2014).
By “carrying information” we means that, for example, we should be able to
extract information about richness or evenness from our index. One way of
extracting such information is to decompose the index additively or
multiplicatively into components that can be interpreted in a biologically
meaningful way; see for example discussions of $\alpha$\- $\beta$\- and
$\gamma$\- diversity (Jost, 2007; Anderson et al., 2011). A priori it is not
clear why such a decomposition would exist, whether it has to be unique, and
in cases of non-uniqueness, what are the conditions for optimality of a
decomposition.

As an example of this approach, van Dam (van Dam, 2019) has recently proposed
a straightforward decomposition of the Leinster–Cobbold (LC) index (Leinster
and Cobbold, 2012). In a sense, our work below is a generalisation of the work
of van Dam, which uses intrinsic properties of the LC index to remove an
important bias in van Dam’s decomposition with intriguing and far-reaching
consequences.

The structure of the paper is as follows. In Section 2 we discuss the
definitions of richness, evenness and similarity using (Chao et al., 2014;
Chiu et al., 2014; Daly et al., 2018; Gregorius and Gillet, 2022). In Section
3 we collect the required information about the LC index following (Leinster
and Cobbold, 2012; Leinster and Meckes, 2016); it subsumes many other
diversity indices such Rao’s index that is widely used in functional ecology
(Rao, 1982; Ricotta and Moretti, 2011). In Section 4 we present van Dam’s and
then our decomposition and its consequences. Finally, in Section 6 we discuss
the relation of our work to that or Chao and Ricotta (Chao and Ricotta, 2019),
extensions and open problems.

##  2\. Diversity components

###  2.1. Notation

First of all, we need to establish notation. Everywhere below we assume that
the number of species in a community is fixed at $n>1$.

We will use ${\boldsymbol{p}}=(p_{1},\,\ldots,p_{n})$ to denote the vector of
relative abundances, and let

(1) $\Delta(n)=\\{{\boldsymbol{x}}\in{\mathbb{R}}^{n},\;|\;x_{k}\geq
0,\;\;k=1,\ldots,n,\;\;\sum_{k=1}^{n}x_{k}=1\\}$

be the standard $n-1$-simplex in ${\mathbb{R}}^{n}$.

######  Remark 1.

It has to be emphasised that admissible relative abundance vectors
${\boldsymbol{p}}$ take values in $\Delta(n)^{\circ}$, the interior of
$\Delta(n)$:

(2)
$\Delta(n)^{\circ}=\\{{\boldsymbol{x}}\in{\mathbb{R}}^{n},\;|\;x_{k}>0,\;\;k=1,\ldots,n,\;\;\sum_{k+1}^{n}x_{k}=1\\},$

which is not a closed set in ${\mathbb{R}}^{n}$; the consequence of that is
that there are converging sequences
$\left({\boldsymbol{p}}_{m}\right)_{m=1}^{\infty}$ in $\Delta(n)^{\circ}$,
whose limit is contained in the boundary
$\partial\Delta(n)=\Delta(n)\backslash\Delta(n)^{\circ}$; such limits by
necessity correspond to communities with fewer than $n$ species.

We will often use the vector

(3)
${\boldsymbol{p}}_{h}=\left(\frac{1}{n},\ldots,\,\frac{1}{n}\right)\in\Delta(n)^{\circ};$

the subscript $h$ stands for “homogeneous”. ${\boldsymbol{p}}_{h}$ is the
relative abundance of a community where each species is represented equally.
We denote by $M(n)\subset\Delta(n)\backslash\Delta(n)^{\circ}$ the set of
$n$-vectors having one component equal to 1 and the rest equal to zero. Thus,
${\boldsymbol{m}}\in M(n)$ is the relative abundance vector of a monomorphic
community. $\bf{1}$ will stand below for an $n$-vector with all components
equal to $1$.

Next, we need to discuss sets of $n\times n$ matrices. First of all, we will
denote the $n\times n$ identity matrix by $I_{n}$. We will use the notation
$J_{n}$ for the $n\times n$ matrix of ones.

In the present paper, for simplicity, we will be working with ultrametric
matrices; this choice is motivated by the fact that similarity matrices (see
subsection 3.1) constructed using taxonomic trees are necessarily ultrametric
and since using them simplifies the theory of (Leinster and Meckes, 2016). For
more information on ultrametric matrices, please see (Dellacheria et al.,
2014, Ch. 3) and (Leinster, 2013; Leinster and Meckes, 2016). We will denote
that set of all ultrametric $n\times n$ matrices by $\mathcal{U}(n)$ and its
interior by $\mathcal{U}(n)^{\circ}$.

######  Definition 2 (Defn. 3.2 of (Dellacheria et al., 2014)).

A symmetric $n\times n$ matrix $A$ is ultrametric if $A_{i,i}\geq\max_{k\neq
i}A_{ik}$, $i,k\in\\{1,\ldots,n\\}$ and $A_{ik}\geq\min\\{A_{ij},A_{jk}\\}$,
$i,j,k\in\\{1,\ldots,n\\}$.

######  Remark 3.

Note that in (Leinster and Meckes, 2016, Example 12) Leinster and Meckes take
the matrices they call ultrametric to be strictly diagonally dominant. That
would preclude the set of ultrametric matrices from being closed; so in our
definition $J_{n}\in\mathcal{U}(n)$.

###  2.2. Evenness

The concept of evenness (for which see, e.g. (Chao et al., 2014; Chao and
Ricotta, 2019; Gregorius and Gillet, 2022)) and references therein), is rather
problematic. First of all, the terminology is badly chosen as it would
immediately seem that the “most even” population of $n$ species is one for
which the vector of relative abundances is the homogeneous vector
${\boldsymbol{p}}_{h}$, i.e. one where every species is equally represented.
Thus, like the case of richness discussed below, the terminology seems to be
precluding discussion. As rightly pointed in (Gregorius and Gillet, 2022),
such a categorical answer to the question of maximal evenness leaves open the
discussion of what would constitute “maximum unevenness” in
$\Delta(n)^{\circ}$. van Dam (van Dam, 2019) uses instead the concept of
“balance”, which seems to us a better term; this is the concept which, after
defining it properly (see (11)) we will be using below.

###  2.3. Richness

Richness is sometimes summarily defined to be the number of species, see e.g.
(Daly et al., 2018). Our approach below allows us to retain this definition
but at a price. Such a definition is open to the same criticism as the notion
of maximal evenness defined by ${\boldsymbol{p}}_{h}$ that we have discussed
above. It is again “species-centric”, and takes into account only the last
level of taxonomic classification. Below, in Section 4 we suggest how to
introduce a defensible new notion of richness that uses taxonomic information.

###  2.4. Similarity

In this section we discuss the construction of taxonomic similarity matrices
$Z$ for a community with $n$ species.

The usual way of constructing similarity matrices $Z$ which are automatically
ultrametric is to assign distances between different levels of a taxonomic
tree. Then the taxonomic distance $d(i,j)$ between two species is the sum of
distance from the nodes corresponding to these species to the first common
node, and then one puts $Z_{ij}=e^{-d(i,j)}$ or if the maximal distance in the
tree has been normalised to 1, one could put $Z_{ij}=1-d(i,j)$. As a example,
consider the tree in Figure 1

$F$ $G_{1}$ $S_{1}$ $S_{2}$ $G_{2}$ $S_{3}$

Figure 1. A taxonomic tree

and set the species-genus and the genus-family distance to be $0.3$. If we use
the additive recipe, we get

(4) $Z_{1}=\left(\begin{array}[]{ccc}1&0.7&0.4\\\ 0.7&1&0.4\\\
0.4&0.4&1\end{array}\right)$

##  3\. The Leinster–Cobbold diversity index

In their influential paper (Leinster and Cobbold, 2012), Leinster and Cobbold
introduced a far-reaching generalisation of Hill numbers, for discussions of
which see (Chiu et al., 2014; Hill, 1973), the LC index. For details on its
properties, see (Leinster and Cobbold, 2012; Leinster and Meckes, 2016); here
we just collect the bare minimum in the framework of taxonomic (ultrametric)
similarity matrices.

###  3.1. Definition of the LC index

As in the definition of Hill numbers, below $q\in[0,\infty)$ is the
sensitivity parameter, measuring the importance given to rare species. Then
for a community of $n$ species with relative abundance vector
${\boldsymbol{p}}$ and (ultrametric) similarity matrix $Z$, we have

######  Definition 4.

The LC diversity of order $q$ is

(5)
$F(Z,{\boldsymbol{p}},q):=\left(\sum_{i=1}^{n}p_{i}\left(Z{\boldsymbol{p}}^{T}\right)_{i}^{q-1}\right)^{1/(1-q)}.$

Note that (Leinster and Cobbold, 2012) use a different notation, similar to
the Hill number notation in the literature; they denote the right-hand side of
(5) by ${}^{q}D^{Z}({\boldsymbol{p}})$. We prefer the notation used here as it
clearly shows functional dependencies and allows easy generalisation, which we
discuss briefly in Section 6. We collect the required properties of the LC
index in the proposition below and in subsection 3.2.

######  Proposition 5.

Let ${\boldsymbol{p}}\in\Delta(n)^{\circ}$. $Z\in U(n)$. Then

  1. [(a)]

  2. (1)

$F({\boldsymbol{p}},Z,q)$ is a monotone decreasing function of $q$;

  3. (2)

$F({\boldsymbol{p}},Z,q)<F({\boldsymbol{p}},I_{n},q)$ for all $q$ if $Z\neq
I_{n}$;

  4. (3)

$F({\boldsymbol{p}}_{h},I_{n},q)=n$ for all $q$;

  5. (4)

$F({\boldsymbol{p}},J_{n},q)=1$ for all $q$.

For proofs of (a) and (b) please see (Leinster and Cobbold, 2012); the rest
are immediate.

Following (Leinster and Meckes, 2016), we now discuss the concept of a
_maximally balanced abundance vector_ for a community of $n$ species with an
ultrametric similarity matrix $Z$.

###  3.2. A crucial property of the LC index

Using only ultrametric taxonomic similarity matrices simplifies the
presentation considerably. For the more general case where the similarity
matrix is simply a symmetric matrix with positive elements, see (Leinster and
Meckes, 2016). The results of that paper have not, in our opinion, been
sufficiently seriously considered by the biodiversity community.

We present two theorems from (Leinster and Meckes, 2016). First of all we have
the following existence and uniqueness result for maximisers of the LC
diversity index.

######  Theorem 6.

For each $Z\in\mathcal{U}(n)^{\circ}$ there exists a unique abundance vector
${\boldsymbol{p}}^{*}\in\Delta(n)^{\circ}$ that maximises
$F(Z,{\boldsymbol{p}},q)$ for every value of $q\in[0,\infty)$.

######  Definition 7.

Given $Z\in\mathcal{U}(n)^{\circ}$, we call the corresponding abundance vector
${\boldsymbol{p}}^{*}$ the maximally balanced abundance vector.

This is the vector that corresponds to ${\boldsymbol{p}}_{h}$ that arises in
theories that do not take into account taxonomic similarity.

Computing the maximally balanced vector in the case of ultrametric similarity
matrices is a simple matter of solving a system of linear equations and
normalising. If the similarity matrix is not ultrametric, the situation is
more complex; see (Leinster and Meckes, 2016) for details.

######  Theorem 8.

Given $Z\in\mathcal{U}(n)^{\circ}$, ${\boldsymbol{p}}^{*}$ is given by

$p^{*}_{i}=\frac{w_{i}}{\sum_{j=1}^{n}w_{j}},$

where ${\boldsymbol{w}}$ solves the system of equations
$Z{\boldsymbol{w}}={\boldsymbol{1}}$, where ${\boldsymbol{1}}$ is a column
vector of ones.

Note that (Leinster and Meckes, 2016, Lemma 6) provides an alternative way of
computing ${\boldsymbol{p}}^{*}$.

######  Definition 9.

A taxonomic tree will be called taxonomically equilibrated if
${\boldsymbol{p}}^{*}={\boldsymbol{p}}_{h}$.

Of course we have

######  Proposition 10.

If at each level of the tree all the nodes have the same degree, the taxonomic
tree is equilibrated.

The converse of Proposition 10 does not hold, i. e. there are taxonomic graphs
that do not satisfy the conditions of that proposition, for which a metric
$d(\cdot,\cdot)$ can be assigned such that the resulting
${\boldsymbol{p}}^{*}$ is ${\boldsymbol{p}}_{h}$. An example is provided by
the following tree:

$O$ $F_{1}$ $G_{1}$ $S_{1}$ $S_{2}$ $G_{2}$ $S_{3}$ $S_{4}$ $F_{2}$ $G_{3}$
$S_{5}$ $S_{6}$ $S_{7}$

Figure 2. A taxonomic tree that allows
${\boldsymbol{p}}^{*}={\boldsymbol{p}}_{h}$.

It is not hard to show that the assignment of species-genus, genus-family and
family-order distances of $0.25$ and using the additive recipe, results in a
similarity matrix for which ${\boldsymbol{p}}^{*}={\boldsymbol{p}}_{h}$. Thus
there is a trichotomy of taxonomic trees: those in Proposition 10 is which
${\boldsymbol{p}}^{*}={\boldsymbol{p}}_{h}$ holds for every assignment of
distances; those where such assignments can be chosen, as in Figure 2, and
such that no assignment of distances results in a homogeneous maximally
balanced abundance vector; an example of such a tree is in Figure 3.

$O$ $F_{1}$ $G_{1}$ $S_{1}$ $G_{2}$ $S_{2}$ $S_{3}$ $F_{2}$ $G_{3}$ $S_{4}$
$S_{5}$

Figure 3. A taxonomic tree for which
${\boldsymbol{p}}^{*}\neq{\boldsymbol{p}}_{h}$ is guaranteed.

We will discuss this trichotomy in more detail in (Chen and Grinfeld, 2023).

######  Remark 11.

Compared to the diversity metrics proposed by Chao et al. (Chao et al., 2014),
the LC index has the flexibility to take into account taxonomic, phylogenetic,
and functional diversity simultaneously. However in that case the resulting
similarity matrices are no longer ultrametric and we leave that more general
case for future study. Also note that the above trichotomy of taxonomic trees
would not necessarily exit if there were a canonical way of constructing
similarity matrices.

##  4\. An unbiased decomposition scheme

We are now ready to propose a decomposition scheme for the LC index. It is
best to start with the decomposition scheme proposed by van Dam (van Dam,
2019) and see why it has to be modified. van Dam writes

(6)
$F({\boldsymbol{p}},Z,q)=\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}},I_{n},q)}\cdot\frac{F({\boldsymbol{p}},I_{n},q)}{F({\boldsymbol{p}}_{h},I_{n},q)}\cdot
F({\boldsymbol{p}}_{h},I_{n},q).$

The first fraction is clearly a measure of dissimilarity, while the second
fraction is a measure of balance; of course
$F({\boldsymbol{p}}_{h},I_{n},q)=n$, so the last term in the right-hand side
is a richness. From many points of view, this is a good decomposition as the
two fractions always lie in the interval $(1/n,1]$. Consult Proposition 5 to
see that the infimum $1/n$ is never reached. The problem here is with the
definition of the measure of balance, as it does not take into account the
taxonomic similarity matrix $Z$ while the dissimilarity measure uses
information from both ${\boldsymbol{p}}$ and $Z$. We call such a decomposition
asymmetrically biased. A different asymmetrically biased decomposition is
given by

(7)
$F({\boldsymbol{p}},Z,q)=\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}}_{h},Z,q)}\cdot\frac{F({\boldsymbol{p}}_{h},Z,q)}{F({\boldsymbol{p}}_{h},I_{n},q)}\cdot
F({\boldsymbol{p}}_{h},I_{n},q).$

In this decomposition the first fraction is a measure of balance, the second a
measure of dissimilarity, while as before, the last term is richness. Of
course here it is the measure of dissimilarity that is asymmetrically biased.
A possibility that might be considered of multiplying (6) and (7) and taking
the square root. That will give us a decomposition which we would call
unbiased (though it could be also called “symmetrically biased”) . However,
there is an additional problem in (7) which is that the first fraction in the
right-hand side can take values larger than one if for example, $Z$ is not a
similarity matrix of a taxonomically equilibrated tree and
${\boldsymbol{p}}={\boldsymbol{p}}^{*}$. We do not pursue this direction as we
do not see any reason for a relative measure not to take values in $[0,1]$.

The price of ensuring normalisation is having to deal with richness in more
detail. Consider instead of (7), the following decomposition:

(8)
$F({\boldsymbol{p}},Z,q)=\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}}^{*},Z,q)}\cdot\frac{F({\boldsymbol{p}}^{*},Z,q)}{F({\boldsymbol{p}}^{*},I_{n},q)}\cdot
F({\boldsymbol{p}}^{*},I_{n},q).$

It is asymmetrically biased as second factor does not involve
${\boldsymbol{p}}$. We will discuss the interpretation of
$F({\boldsymbol{p}}^{*},I_{n},q)$ later.

To obtain an unbiased decomposition, we therefore multiply (6) and (8) and
take a square root. The result is

(9)
$F({\boldsymbol{p}},Z,q)=\sqrt{\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}}^{*},Z,q)}\frac{F({\boldsymbol{p}},I_{n},q)}{F({\boldsymbol{p}}_{h},I_{n},q)}}\cdot\sqrt{\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}},I_{n},q)}\frac{F({\boldsymbol{p}}^{*},Z,q)}{F({\boldsymbol{p}}^{*},I_{n},q)}}\cdot\sqrt{nF({\boldsymbol{p}}^{*},I_{n},q)}.$

The last term in the right-hand side can be rewritten as

(10)
$\sqrt{nF({\boldsymbol{p}}^{*},I_{n},q)}=\sqrt{\frac{F({\boldsymbol{p}}^{*},I_{n},q)}{n}}n:=E(Z,q)n.$

The term $E(z,q)$ expresses the lack of equilibration in the taxonomic tree;
see Definition 9 and Theorem 10.

Putting

(11)
$B({\boldsymbol{p}},Z,q)=\sqrt{\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}}^{*},Z,q)}\frac{F({\boldsymbol{p}},I_{n},q)}{F({\boldsymbol{p}}_{h},I_{n},q)}},$
(12)
$D({\boldsymbol{p}},Z,q)=\sqrt{\frac{F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}},I_{n},q)}\frac{F({\boldsymbol{p}}^{*},Z,q)}{F({\boldsymbol{p}}^{*},I_{n},q)}},$

We finally write our decomposition as

(13)
$F({\boldsymbol{p}},Z,q)=B({\boldsymbol{p}},Z,q)D({\boldsymbol{p}},Z,q)E(Z,q)n.,$

i.e. a product of measures of balance $B({\boldsymbol{p}},Z,q)$, dissimilarity
$D({\boldsymbol{p}},Z,q)$, (lack of) equilibration $E(Z,q)$ and the classical
richness $n$.

Note that by construction both the measure of balance (11) and dissimilarity
(12) are constrained to lie in $[0,1]$ by Proposition 5.

Both the measure of balance and of dissimilarity are geometric means of an
unbiased measure and a biased one. It does not seem possible to find a truly
unbiased decomposition of the LC index, which is the reason we could call the
decomposition (13) symmetrically biased.

Note that though $B({\boldsymbol{p}},Z,q)\geq
1/(\sqrt{F({\boldsymbol{p}}^{*},Z,q)}$ (the right-hand side being independent
of $q$), it is not clear what vector ${\boldsymbol{p}}(q)$ maximises it for a
particular value of $q$. Of course the value $1$ is reached for $q=0$ by the
choice ${\boldsymbol{p}}(0)={\boldsymbol{p}}^{*}$.

$E(Z,q)$ depends on $Z$, as $Z$ defines ${\boldsymbol{p}}^{*}$, and hence
$E(Z,q)$ reflects the structure of the underlying taxonomic tree. Note that in
the case of similarity matrices of taxonomically equilibrated trees, for which
we have ${\boldsymbol{p}}^{*}={\boldsymbol{p}}_{h}$ and hence
$F({\boldsymbol{p}}^{*},I_{n},q)=n$, so that $E(Z,q)=1$, If $Z$ does not
correspond to a taxonomically equilibrated taxonomic tree,
$F({\boldsymbol{p}}^{*},I_{n},q)$ is dependent on $q$.

######  Remark 12.

We could have defined a notion of “richness” by $R(z,q):=E(Z,q)n$, but the
decomposition (13) seems to us more insightful and does not necessitate an
advocacy of a new notion of richness. Definition 9 singles out a class of
taxonomic trees for which the two notions of richness coincide.

##  5\. Practical examples

To illustrate how our decomposition approach differs from the “ABC” approach
suggested by van Dam (van Dam, 2019), we use a simple example in Leinster and
Cobbold (Leinster and Cobbold, 2012) (their Example 3; original data from
(DeVries et al., 1997)). This example has the nice feature that the difference
between the two communities changes the sign when $q$ increases from 0 to 2.

Figure 4. Comparing our approach (“This study”) and van Dam (van Dam, 2019) to
decompose (A) the LC diversity index into (B) Lack of equilibration, (C)
Balance (or Evenness), and (D) Dissimilarity for the example of Charaxinae in
Leinster and Cobbold (Leinster and Cobbold, 2012).

Compared to van Dam (van Dam, 2019), our new decomposition approach yields
different interpretations regarding what aspects of diversity lead to the
differences in diversity estimates between the two communities (Canopy vs.
Understory). When we put more emphasis on rare species ($q<1$), the Canopy
community is more diverse because it has a larger number of species richness
(6 vs. 5) and its species are slightly more dissimilar with each other (Fig.
4A,D). By contrast, when we focus more on abundant species ($q>1$), the
Understory community becomes more diverse because of its greater balance
(evenness) of dominant species (Fig. 4A,C). The difference between our
decomposition approach and van Dam (van Dam, 2019) is that our approach
predicts a much larger difference in Balance than van Dam (van Dam, 2019) (and
a smaller difference in species dissimilarity) when we increasingly focus on
dominant species (larger $q$). Our approach also shows that as $q$ increases,
the Understory community shows a stronger lack of equilibration (i.e.,
deviation of ${\boldsymbol{p}}$ from ${\boldsymbol{p}}^{*}$) than the Canopy
community, albeit this difference is small. In summary, we would interpret
that the greater diversity of the Understory community when we focus on
dominant species is because its dominant species are more balanced than those
in the Canopy community. On the contrary, the interpretation would be that the
dominant species are more dissimilar to each other in the Understory community
if we use the approach of van Dam (van Dam, 2019).

##  6\. Discussion

The LC index uses three “information streams”: the number of species
(richness) $n$, the relative abundance vector ${\boldsymbol{p}}$ and the
similarity matrix $Z$. We could in theory consider a diversity index
$F(c_{1},\;\ldots,c_{m};\,q)$, where $c_{1},\ldots,\;c_{m}$ are information
streams, sources of information about the structure of the community,
expressed as some mathematical objects (vectors, matrices higher order
tensors). We could then follow the decomposition process of Section 4: find
$m!$ biased decompositions, multiply them together and take the $m!$-root.
However, this is already unwieldy in the case of $m=3$. But note that this
procedure is unnecessary as the LC theory has a lot of built-in flexibility in
the definition of similarity matrices. As explained in subsection 2.4, one can
define a similarity matrix by setting $Z_{ij}=e^{-d(i,j)}$, where $d(i,j)$ is
some suitably defined distance between species $i$ and $j$. Hence
incorporating more information streams can be thought about as changing the
distance function $d(\cdot,\cdot)$; in the process of incorporating such
information, such as functional similarity, the ultrametricty of the
similarity matrix is lost; it is possible that the resulting function
$d(\cdot,\cdot)$ will no longer be a metric, becoming more generally a
divergence measure. The point is ${\boldsymbol{p}}$ and (a suitably redefined)
$Z$ dependence of a diversity index is sufficient to incorporate all relevant
information.

In (Chao and Ricotta, 2019), Chao and Ricotta show how to quantify evenness
using divergence measures. It is useful therefore to consider the LC diversity
index (4) in this context. As now we deal with balance as opposed to evenness,
we will denote the resulting balance index by $B$. First of all, let us note
that (11) cannot give rise to an divergence measure-based index of balance as
there is no well defined upper bound for it for all $q$. It is still of course
useful in providing an estimate of the balance contribution to the LC
diversity index. These two statements are not in contradiction.

Clearly, the LC index itself provides a divergence measure-based estimate
(index) of balance via

(14)
${\mathbb{B}}=\frac{F({\boldsymbol{p}},Z,q)-1}{F({\boldsymbol{p}}^{*},Z,q)-1},$

where one could alternatively write $1=F({\boldsymbol{m}},Z,q)$ where
${\boldsymbol{m}}$ is any vector in $M(n)$.

Concerning similarity indices, it again does not seem to be possible to
utilise $D({\boldsymbol{p}},Z,q)$ of (12) to this end. On the other had, the
LC index itself provides a divergence measure based similarity index by

(15)
${\mathbb{S}}=\frac{F({\boldsymbol{p}},I_{n},q)-F({\boldsymbol{p}},Z,q)}{F({\boldsymbol{p}},I_{n},q)-F({\boldsymbol{p}},J_{n},q)}.$

Again, here the value 1 is never reached over $\mathcal{U}(n)^{\circ}$. It is
not hard to show the following proposition:

######  Proposition 13.

The indices ${\mathbb{B}}$, ${\mathbb{S}}$ satisfy all the requirements in
(Chao and Ricotta, 2019).

To conclude, we have proposed a novel decomposition of the LC index. Compared
to a previous version of decomposition due to van Dam (van Dam, 2019), our
approach estimates the balance and dissimilarity of the community more
comprehensively (e.g., we not only estimate dissimilarity for a homogeneous
community but also consider the present vector of relative abundance). As
such, we believe that our inference is more robust when comparing balance and
dissimilarity among communities. In addition, we had by necessity to introduce
a notion of taxonomic tree equilibration (which turns out to be an important
concept in our on-going work on quantifying un-evenness (Chen and Grinfeld,
2023)), which is another descriptor of a biological community. We advocate the
use of our decomposition (13) as a “maximally unbiased” estimate of
contributions of balance and (dis)similarity to diversity.

## References

  * Anderson et al. (2011) Anderson, M. J., Crist, T. O., Chase, J. M., Vellend, M., Inouye, B. D., Freestone, A. L., Sanders, N., Cornell, H., Comita, L., Davies, K., and Harrison, S. (2011).  Navigating the multiple meanings of $\beta$ diversity: a roadmap for the practicing ecologist.  Ecology Letters, 14:19–28.
  * Chao et al. (2014) Chao, A., Chiu, C.-H., and Jost, L. (2014).  Unifying species diversity, phylogenetic diversity, functional diversity, and related similarity differentiation measures through Hill numbers.  Annual Review of Ecology, Evolution, and Systematics, 45:297–324.
  * Chao and Ricotta (2019) Chao, A. and Ricotta, C. (2019).  Quantifying evenness and linking it to diversity, beta diversity and similarity.  Ecology, 100:e02852.
  * Chen and Grinfeld (2023) Chen, B. and Grinfeld, M. (2023).  An evennesss index based on the Leinster–Cobbold diversity index.  in progress.
  * Chiu et al. (2014) Chiu, C.-H., Jost, L., and Chao, A. (2014).  Phylogenetic beta diversity, similarity, and differentiation measures based on Hill numbers.  Ecological Monographs, 81:21–44.
  * Daly et al. (2018) Daly, A. J., Baetens, J. M., and De Baets, B. (2018).  Ecological diversity: measuring the unmeasurable.  Mathematics, 6:119.
  * Dellacheria et al. (2014) Dellacheria, C., Martinez, S., and San Martin, J. (2014).  Inverse M-matrices and Ultrametric Matrices.  Springer.
  * DeVries et al. (1997) DeVries, P. J., Murray, D., and Lande, R. (1997).  Species diversity in vertical, horizontal, and temporal dimensions of a fruit-feeding butterfly community in an ecuadorian rainforest.  Biological Journal of the Linnean Society, 62(3):343–364.
  * Gregorius and Gillet (2022) Gregorius, H.-R. and Gillet, E. M. (2022).  The concept of evenness/unevenness: less evenness or more unevenness?  Acta Biotheoretica, 70:1–28.
  * Hill (1973) Hill, M. (1973).  Diversity and evenness: a unifying notation and its consequences.  Ecology, 54:427–432.
  * Hillebrand et al. (2018) Hillebrand, H., Blasius, B., Borer, E. T., Chase, J. M., Downing, J. A., Eriksson, B. K., Filstrup, C. T., Harpole, W., Hodapp, D., Larsen, S., and Lewandowska, A. M. (2018).  Biodiversity change is uncoupled from species richness trends: Consequences for conservation and monitoring.  Journal of Applied Ecology, 55(1):169–184.
  * Jost (2007) Jost, L. (2007).  Partitioning diversity into independent alpha and beta components.  Ecology, 88:2427–2439.
  * Leinster (2013) Leinster, T. (2013).  The magnitude of metric spaces.  Documenta Mathematica, 18:857–905.
  * Leinster and Cobbold (2012) Leinster, T. and Cobbold, C. (2012).  Measuring diversity: the importance of species similarity.  Ecology, 93:477–489.
  * Leinster and Meckes (2016) Leinster, T. and Meckes, M. W. (2016).  Maximizing diversity in biology and beyond.  Entropy, 18:88.
  * Maureaud et al. (2019) Maureaud, A., Hodapp, D., van Denderen, P. D., Hillebrand, H., Gislason, H., Spaanheden Dencker, T., Beukhof, E., and Lindegren, M. (2019).  Biodiversity–ecosystem functioning relationships in fish communities: biomass is related to evenness and the environment, not to species richness.  Proceedings of the Royal Society B, 286(1906):20191189.
  * Purvis and Hector (2000) Purvis, A. and Hector, A. (2000).  Getting the measure of biodiversity.  Nature, 405(6783):212–219.
  * Rao (1982) Rao, C. R. (1982).  Diversity and dissimilarity coefficients: a unified approach.  Theoretical population biology, 21(1):24–43.
  * Ricotta and Moretti (2011) Ricotta, C. and Moretti, M. (2011).  CWM and Rao’s quadratic diversity: a unified framework for functional ecology.  Oecologia, 167(1):181–188.
  * Stuart-Smith et al. (2013) Stuart-Smith, R. D., Bates, A. E., Lefcheck, J. S., Duffy, J. E., Baker, S. C., Thomson, R. J., Stuart-Smith, J. F., Hill, N. A., Kininmonth, S. J., Airoldi, L., and Becerro, M. A. (2013).  Integrating abundance and functional traits reveals new global hotspots of fish diversity.  Nature, 501(7468):539–542.
  * van Dam (2019) van Dam, A. (2019).  Diversity and its decomposition into variety, balance and disparity.  Royal Society open science, 6:90452.