# Morse index and determinant of block Jacobi matrices via optimal control

Stefano Baranzini 111Università degli studi di Torino, Via Carlo Alberto 10,
Torino, Italy  Ivan Beschastnyi 222Univesidade de Aveiro, Campus Santiago,
Aveiro, Portugal

###### Abstract

We describe the relation between block Jacobi matrices and minimization
problems for discrete time optimal control problems. Using techniques
developed for the continuous case, we provide new algorithms to compute
spectral invariants of block Jacobi matrices. Some examples and applications
are presented.

##  1 Introduction

###  1.1 Motivation

The goal of this paper is to explore an interesting connection between block
Jacobi matrices and a class of discrete optimal control problems. This gives
effective algorithms for computing the determinant, the negative inertia index
and more generally the number of eigenvalues smaller than some threshold
$\lambda^{*}\in\mathbb{R}$ of large block Jacobi matrices.

Recall that a block Jacobi matrix $\mathcal{I}$ is a matrix of the form

${\mathcal{I}}=\begin{pmatrix}S_{1}&R_{1}&0&\dots&0&0\\\
R_{1}^{t}&S_{2}&R_{2}&\dots&0&0\\\ 0&R_{2}^{t}&S_{3}&\dots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\dots&S_{N-1}&R_{N-1}\\\
0&0&0&\dots&R^{t}_{N-1}&S_{N}\end{pmatrix},$ (1)

where $N\in\mathbb{N}$, $S_{k}$ are symmetric matrices of order $n$. Usually
it is also assumed that $R_{k}$ are positive symmetric matrices of the same
order. Here, we only require them to be non-degenerate, i.e. $R_{k}\in
GL(n,\mathbb{R})$.

Jacobi matrices find applications in numerical analysis [11], statistical
physics [12], knot theory [8] and many other areas. Moreover, any symmetric
matrix can be put in a tridiagonal form using the Householder method [7].
Therefore, understanding their spectral properties is an extremely important
topic.

In this article we establish a correspondence between Jacobi matrices and
discrete linear quadratic regulator (LQR) problem. Then we show how the
general machinery of optimal control theory allows us to compute the index and
the determinant of such matrices. Similar formulas for determinants also
appeared in [15], see also [16, 17] for a geometric approach to the problem.

The LQR approach nevertheless has its benefits. In particular, as a small
application of our method we give a novel proof of the generalized Euler
identity

$\frac{a\sin b}{b\sinh
a}=\prod_{j=1}^{\infty}\left(1-\frac{a^{2}+b^{2}}{a^{2}+(\pi j)^{2}}\right)$
(2)

using only methods of linear algebra and basic analysis.

In order to state the main results, let us introduce the necessary notations.
Consider the following discrete differential equation:

$x_{k+1}-x_{k}=A_{k}x_{k}+B_{k}u_{k},\qquad k=0,\dots,N+1,$ (3)

Where $A_{k}$ and $B_{k}$ are $n\times n$ matrix, $B_{k}$ are invertible and
$x_{k},u_{k}\in\mathbb{R}^{n}$. Variables $x_{k}$ are called _state_ variables
while $u_{k}$ are the _controls_. We denote by $x=(x_{k})$ and $u=(u_{k})$ the
vectors having $x_{k}$ and $u_{k}$ respectively as $k$-th component. Note that
we will often consider them as elements of the space of functions on a finite
set of points with values in $\mathbb{R}^{n}$, clearly isomorphic to
$\mathbb{R}^{n(N+2)}$ and $\mathbb{R}^{n(N+1)}$ correspondingly. On the space
of solutions of eq. 3 with fixed initial condition $x_{0}$ we consider the
following quadratic functional given by:

$\mathcal{J}(u)=\frac{1}{2}\sum_{i=0}^{N}\left(|u_{i}|^{2}-\langle
Q_{i}x_{i},x_{i}\rangle\right),$ (4)

where $Q_{i}$ are some symmetric $n\times n$ matrices.

The problem of minimizing $\mathcal{J}$ over all possible discrete
trajectories of (3) is the classical linear quadratic regulator problem, which
is one of the central problems in optimal control theory. Since we assume that
$B_{k}$ are invertible, we can resolve the constraints (3) with respect to
$u_{k}$ and plug them into (4). This yields

$\mathcal{J}(x)=\frac{1}{2}\left(\sum_{i=0}^{N}|B_{i}^{-1}(x_{i+1}-(1+A_{i})x_{i})|^{2}-\sum_{i=0}^{N}\langle
Q_{i}x_{i},x_{i}\rangle\right).$

If we assume additionally that $x_{0}=0=x_{N+1}$ and write down the matrix
associated to the quadratic form above, we find that it is given by a block
Jacobi matrix (1). Indeed, set

$\Gamma_{k}^{-1}:=(B_{k}^{-1})^{t}B_{k}^{-1},$

which is a symmetric non degenerate matrix. Then the functional $2\mathcal{J}$
is the represented by the symmetric matrix (1) with:

$\displaystyle S_{k}$
$\displaystyle=\Gamma_{k-1}^{-1}+(1+A_{k})^{t}\Gamma_{k}^{-1}(1+A_{k})-Q_{k},$
$\displaystyle R_{k}$ $\displaystyle=-(1+A_{k})^{t}\Gamma_{k}^{-1},$ (5)
$\displaystyle 1\leq k\leq N-1.$

We now see that to each LQR problem (3)-(4) we can associate a Jacobi matrix.
The converse is also true. However, there are several LQR problems which
correspond to the same Jacobi matrix. We can get rid of this ambiguity by
choosing fixed invertible $B_{k}$, $0\leq k\leq n$ and an arbitrary matrix
$A_{0}$. Thus to every LQR problem (3)-(4) one can associate a Jacobi matrix
and vice versa. Fixing $B_{k}$, $0\leq k\leq n$ and $A_{0}$ beforehand makes
this correspondence an affine bijection.

###  1.2 Main results

Through out the paper we will make the following assumption. It is equivalent
to requiring $R_{k}\in GL(n,\mathbb{R})$ in (5).

######  Assumption 1.

The matrices $B_{k}$ and $1+A_{k}$ in (3) are invertible for every $k$.

Similarly to the case of continuous optimal control problems, one can
characterize the extremal points of the functional (4) under the constraints
(3) as solutions to a system of difference equations. To be more precise
define the following matrices for $0\leq i\leq j\leq N+1$

$P_{i}^{j}:=\prod_{r=i}^{j-1}(1+A_{r}),\quad P_{j}^{j}=1.$

and for $0\leq i<j\leq N+1$ we take:

$P_{j}^{i}=(P_{i}^{j})^{-1}.$

Next, we define $2n\times 2n$ matrices:

$\displaystyle M_{k}=\begin{pmatrix}(P^{k}_{k+1})^{t}&0\\\
\Gamma_{k}(P^{k}_{k+1})^{t}&P^{k+1}_{k}\end{pmatrix}\begin{pmatrix}1&-Q_{k}\\\
0&1\end{pmatrix}$

and

$\Phi_{k}=\begin{pmatrix}\Phi^{1}_{k}&\Phi^{2}_{k}\\\
\Phi^{3}_{k}&\Phi^{4}_{k}\end{pmatrix}:=\prod_{i=0}^{k-1}M_{i}.$

Matrices $\Phi_{k}$ correspond to the flow of the extremal equations that we
will derive in Section 2 using the Lagrange multiplier rule. Now we can
formulate the two main results of this paper.

######  Theorem 1.

Let $\mathcal{J}$ be as in (4) and consider constraints as in (3). The
negative inertia index of $\mathcal{J}$ restricted to solutions of (3) with
$x_{0}=x_{N+1}=0$ is given by the following formula:

$ind^{-}\mathcal{J}=\sum_{k=1}^{N}ind^{-}\mathcal{Q}_{k}^{\perp}+\dim\ker\Phi^{3}_{k},$
(6)

where

$\mathcal{Q}_{k}^{\perp}=(\Phi_{k}^{3})^{t}\left((P^{k+1}_{k})^{t}\Gamma_{k}^{-1}(P^{k+1}_{k})-Q_{k}\right)\Phi_{k}^{3}+(\Phi_{k}^{3})^{t}\Phi_{k}^{1}$

Notice that we have reduced the problem of computing the negative inertia
index of a $nN\times nN$ matrix to computing the inertia of $N$ matrices of
dimension $n\times n$, which greatly reduces the dimensionality of the
problem. In the same spirit, we also prove a formula for the determinant.

######  Theorem 2.

Let $\mathcal{J}$ be as in (4) and consider constraints as in (3). The
determinant of $\mathcal{J}$ restricted to solutions of (3) with
$x_{0}=x_{N+1}=0$ is given by the following formula:

$\det(\mathcal{J})=2^{-nN}\det(\Phi_{3})\det(P_{0}^{N+1})\prod_{k=0}^{N}\det\Gamma_{k}^{-1}.$
(7)

Let us rephrase the two statements in terms of Jacobi matrices of the form
(1). We can choose $A_{0}=0$ and $\Gamma_{k}=1$ for all $0\leq k\leq N$ and
obtain the following relations from eq. 5 for $1\leq k\leq N-1$:

$\begin{cases}Q_{k}=1-S_{k}+R_{k}R^{t}_{k},\\\
P_{k}^{k+1}=-R_{k}^{t},\end{cases}$
$M_{k}=\begin{pmatrix}-R^{-1}_{k}&R_{k}^{-1}(1-S_{k}+R_{k}R_{k}^{t})\\\
-R_{k}^{-1}&R_{k}^{-1}(1-S_{k})\end{pmatrix}.$

######  Corollary 1.

Index and determinant of the matrix $\mathcal{I}$ defined in eq. 1 can be
computed as:

$\displaystyle ind^{-}\mathcal{I}$
$\displaystyle=\sum_{k=1}^{N}ind^{-}\mathcal{Q}_{k}^{\perp}+\dim\ker\Phi^{3}_{k},$
$\displaystyle\det(\mathcal{I})$
$\displaystyle=(-1)^{nN}\det(\Phi_{3})\prod_{k=1}^{N}\det R_{k}.$

where
$\mathcal{Q}_{k}^{\perp}=(\Phi_{k}^{3})^{t}\left(S_{k}-1\right)\Phi_{k}^{3}+(\Phi_{k}^{3})^{t}\Phi_{k}^{1}$.

######  Remark 1.

Clearly, applying formula (6) with $\tilde{Q}_{k}=Q_{k}+\lambda^{*}$ instead
of $Q_{k}$ (or equivalently $\tilde{S}_{k}=S_{k}-\lambda^{*}$ instead of
$S_{k}$), gives a formula to compute the number of eigenvalues of
$\mathcal{J}$ (respectively of $\mathcal{I}$) smaller than
$\lambda^{*}\in\mathbb{R}$. For some oscillation results in this setting, see
[5].

The techniques used in our proofs are deeply connected to symplectic geometry.
One can extend them to much more general settings [3, 4, 6]. We have chosen to
keep the exposition as simple as possible and minimize the use of symplectic
language in this work. In this way, both formulas (6) and (7) can be applied
directly.

The article has the following structure. In Section 2 we derive the
Hamiltonian system for the extremal curves using the end-point map and
Lagrange multiplier rule. In Section 3 we prove Theorem 1 and in Section 4
Theorem 2. Finally, in Section 5, we consider some examples and prove the
generalized Euler identity (2).

##  2 Eigenvalues of Jacobi matrices and the Lagrange multiplier rule

In order to compute the index and the determinant of a Jacobi matrix
${\mathcal{I}}$ in (1) we will use the flow of a certain system of difference
equations. First, though, we need some notions from optimal control theory.

The first object we introduce is the endpoint map:

$E^{k+1}:u\mapsto x_{k+1}(u),\qquad 0\leq k\leq N,$

This map takes a control $u$ and gives the solution of the differential
equation (3) with $x_{0}=0$ at step ${k+1}$. We can write this map explicitly
by iterating (3).

$E^{k+1}(u)=x_{k+1}(u_{0},\dots,u_{k})=\sum_{j=1}^{k+1}P_{j}^{k+1}B_{j-1}u_{j-1}.$
(8)

In particular, from this formula, it is clear that the value of $E^{k+1}(u)$
depends only on the first $k$ components of the control $u$. It makes sense,
thus, to introduce the following filtration (i.e. _flag_ of subspaces) in the
space of controls. Set $U:=\mathbb{R}^{nN}$ the space of controls and define:

$U_{k}:=\\{(u_{0},u_{1},\dots,u_{k},0,\dots,0):u_{j}\in\mathbb{R}^{n}\\}.$ (9)

Denote by $pr_{l}$ the orthogonal projection on $U_{l}$ with respect to the
standard Euclidean scalar product. We have:

$E_{x_{0}}^{k+1}(u)=E_{x_{0}}^{k+1}(pr_{l}u)\quad\forall\,l\geq k.$

Under 1, we can resolve the constraints $u\in\ker E^{k+1}\cap U_{k}$ quite
easily. This will play an important role in all of the proofs. Indeed,
consider the map

$F_{k}:U_{k-1}\to\mathbb{R}^{n},$

defined as

$F_{k}(u):=-B_{k}^{-1}\sum_{j=1}^{k}P^{k+1}_{j}B_{j-1}u_{j-1}.$ (10)

It is straightforward to check, using eq. 8, that:

$(u_{0},\dots,u_{k})\in\ker E^{k+1}\cap U_{k}\iff u_{k}=F_{k}(u).$ (11)

We can identify $F_{k}$ with a row block-matrix of size $n\times nk$.

As we have seen in the previous section, a Jacobi matrix can be interpreted as
the functional $\mathcal{J}$ in eq. 4 restricted to $\ker E^{N+1}$. We would
like to find the eigenvalues of this restriction, i.e. of the corresponding
Jacobi matrix ${\mathcal{I}}$.

By definition $\tau\in\mathbb{R}$ is an eigenvalue of $\mathcal{J}$ if and
only if the bilinear form $\mathcal{J}-\tau|u|^{2}$ is degenerate on $\ker
E^{N+1}$. Since the form $\mathcal{J}(u)-\tau|u|^{2}$ is quadratic and $\ker
E^{N+1}$ is a linear space, this is equivalent to finding critical points of
$\mathcal{J}(u)-\tau|u|^{2}$ on $\ker E^{N+1}$. This is a typical task in
constrained optimization and can be done via Lagrange multiplier rule. Notice
that $\tau=\frac{1}{2}$ is not an eigenvalue of $\mathcal{J}$ if and only if
all $Q_{i}$, $1\leq i\leq N$ are non degenerate.

Let us multiply $\mathcal{J}-\tau$ by $(1-2\tau)^{-1}$, operation that clearly
preserves the kernel. After we introduce a new parameter $s=1/(1-2\tau)$ we
obtain a family of functionals

$\displaystyle\mathcal{J}_{s}(u):$
$\displaystyle=s\mathcal{J}(u)-\frac{s-1}{2}|u|^{2}$
$\displaystyle=\frac{1}{2}\left(\sum_{i=0}^{N}|u_{i}|^{2}-s\sum_{i=1}^{N}\langle
x_{i},Q_{i}x_{i}\rangle\right).$

The following result is a discrete version of the Pontryagin maximum principle
for LQR problems, which characterizes critical points of the constrained
variational problem as solutions to a Hamiltonian system.

######  Lemma 1.

The quadratic form $\mathcal{J}_{s}$ is degenerate on $\ker E^{N+1}$ if and
only if the boundary value problem

$\begin{cases}x_{k+1}=P^{k+1}_{k}x_{k}+B_{k}B_{k}^{t}\lambda_{k+1},\\\
\lambda_{k}=(P^{k+1}_{k})^{t}\lambda_{k+1}+sQ_{k}x_{k},\end{cases}$ (12)

with $x_{0}=0$, $x_{N+1}=0$ has a non-trivial solution. Moreover the dimension
of the space of solutions of eq. 12 and of $\ker\mathcal{J}_{s}|_{\ker
E^{N+1}}$ coincide.

The former quadratic form is degenerate if and only if $(1-s)/s$ is an
eigenvalue of the restriction to $\ker dE^{N+1}$ of:

$\sum_{i=0}^{N}|v_{i}|^{2}-\sum_{j=1}^{N+1}(E^{i}(v))^{t}Q_{i}E^{i}(v).$

###### Proof.

If $u$ is critical point of $\mathcal{J}_{s}$ restricted to $\ker E^{N+1}$,
there exists $\lambda_{N+1}\in\mathbb{R}^{n}$ such that

$d_{u}\mathcal{J}^{N}=\lambda_{N+1}^{t}d_{u}E^{N+1}.$ (13)

Since $E^{N+1}$ is linear, we have $d_{u}E^{N+1}=E^{N+1}$. We will often use
the latter notation in the formulas that follow.

In order to obtain a discrete analogue of a Hamiltonian system of PMP, we need
to introduce for $1\leq k\leq N$ a family of functionals

$\mathcal{J}_{s}^{k}(u):=\frac{1}{2}\left(\sum_{i=0}^{k}|u_{i}|^{2}-s\sum_{i=1}^{k}\langle
x_{i},Q_{i}x_{i}\rangle\right)$

and their differentials:

$d_{u}\mathcal{J}^{k}_{s}(v)=\sum_{i=0}^{k}\langle
u_{i},v_{i}\rangle-s\sum_{i=1}^{k}x_{i}^{t}Q_{i}E^{i}(v).$ (14)

Functionals $\mathcal{J}^{k}_{s}$ differ from $\mathcal{J}_{s}$ only by the
range of summation and clearly $\mathcal{J}^{N}_{s}=\mathcal{J}_{s}$. We look
for $\lambda_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N+1$ which satisfy:

$d_{u}\mathcal{J}^{k}_{s}=\lambda_{k+1}^{t}d_{u}E^{k+1},$ (15)

Using the explicit formula (8) for the end-point map we find the recurrence
relation

$E^{k+1}(v)=P^{k+1}_{k}E^{k}(pr_{k-1}v)+B_{k}v_{k}$ (16)

for all $v\in U_{k}$. From here, we can already derive the expression for the
$k-th$ component of the control. Indeed, assume that $v_{i}=0$ for all $i<k$.
In this case

$E^{i}(pr_{i-1}v)=0,\quad\forall i\leq k.$

Hence, if we substite such a control $v$ in (15), then using (14) and (16) we
obtain the control law

$u_{k}=B_{k}^{t}\lambda_{k+1},$ (17)

which we can plug in (3).

The next step is to obtain a discrete equation for $\lambda_{k}$. To do this
we compare the multipliers $\lambda_{k}$ and $\lambda_{k+1}$ by subtracting
(15) from the same formula with $k$ replaced by $k-1$. Using $v\in
U_{k-1}\subset U_{k}$ as a test variation and formulas (14) and (16) we find
that:

$\displaystyle\lambda_{k+1}^{t}E^{k+1}(v)-\lambda_{k}^{t}E^{k}(v)=d_{u}\mathcal{J}^{k}_{s}(v)-d_{u}\mathcal{J}^{k-1}_{s}(v)$
$\displaystyle\iff\left(\lambda^{t}_{k+1}P^{k+1}_{k}-\lambda^{t}_{k}+sx_{k}^{t}Q_{k}\right)E^{k}(v)=0.$

Since by our assumption $B_{k}$ are invertible for all $0\leq k\leq N$, the
end-point map $E^{k}$ is surjective, as can be seen from the explicit
expression (8). Thus the term in the bracket must vanish.

Collecting everything gives eq. 12. The boundary conditions for this system
come form the fact that we have to add the equation $E^{N+1}(u)=0=x_{N+1}$ and
that we are assuming $x_{0}=0$.

Now, notice that the initial covector $\lambda_{0}$ of the lift uniquely
determines the whole trajectory and the control. This means that the
correspondence $\lambda_{0}\mapsto u$ is a bijection between the space of
solutions of the boundary value problem and the kernel of $\mathcal{J}_{s}$. ∎

Under 1, we can rewrite system (12) as a forward equation. To do so, recall
that $\Gamma_{k}=B_{k}B_{k}^{t}$, multiply the second equation by
$(P^{k}_{k+1})^{t}$ and plug in the new expression for $\lambda_{k+1}$ into
the first equation. This gives

$\begin{cases}x_{k+1}=(P^{k+1}_{k}-s\Gamma_{k}(P^{k}_{k+1})^{t}Q_{k})x_{k}+\Gamma_{k}(P^{k}_{k+1})^{t}\lambda_{k},\\\
\lambda_{k+1}=(P^{k}_{k+1})^{t}\lambda_{k}-s(P^{k}_{k+1})^{t}Q_{k}x_{k}.\end{cases}$
(18)

We thus have

$\begin{pmatrix}\lambda_{k+1}\\\
x_{k+1}\end{pmatrix}=M_{k}(s)\begin{pmatrix}\lambda_{k}\\\
x_{k}\end{pmatrix},$

where

$\displaystyle M_{k}(s)=\begin{pmatrix}(P^{k}_{k+1})^{t}&0\\\
\Gamma_{k}(P^{k}_{k+1})^{t}&P^{k+1}_{k}\end{pmatrix}\begin{pmatrix}1&-sQ_{k}\\\
0&1\end{pmatrix}.$

Both matrices in the product are symplectic, which makes $M_{k}(s)$ symplectic
as well.

Define the flow up to the point $k$ as

$\Phi_{k}(s)=\begin{pmatrix}\Phi^{1}_{k}(s)&\Phi^{2}_{k}(s)\\\
\Phi^{3}_{k}(s)&\Phi^{4}_{k}(s)\end{pmatrix}:=\prod_{i=0}^{k-1}M_{i}(s).$

When $k=N+1$, we just write $\Phi(s)=\Phi_{N+1}(s)$. Also in accordance with
the notations given in the introduction, we write for $s=1$:

$\mathcal{J}^{k}:=\mathcal{J}^{k}_{1},\qquad\Phi_{k}:=\Phi_{k}(1),\qquad
M_{k}:=M_{k}(1).$

######  Remark 2.

We can reformulate Lemma 1 and the boundary value problem (12) in terms of
$\Phi_{k}(s)$. We are looking for solutions of (12) with boundary values
$x_{0}=x_{k+1}=0$. Writing as above $\Phi_{k+1}(s)$ as a block matrix shows
that the existence of a solution is equivalent to the vanishing of the
determinant of the block $\Phi^{3}_{k+1}$:

$\ker\mathcal{J}^{k}_{s}|_{\ker E^{k+1}}\neq 0\iff\det\Phi^{3}_{k+1}(s)=0.$

##  3 A recursive formula for the index

In this section we prove Theorem 1. The tools we will employ are essentially
taken from linear algebra. The same kind of ideas has been applied in [4] to
reformulate problems of intersection theory for paths in the Lagrange
Grassmannian in terms of finite dimensional linear algebra. The classical
approach to second order minimality conditions rephrases the problem of
counting the number of negative eigenvalues of the Second Variation (the
counter part of our quadratic form $\mathcal{J}$) as counting the intersection
number of a curve in a finite-dimensional manifolds, known as the Lagrangian
Grassmannian, with an appropriate submanifold. For further references on the
continuum approach see for instance [9, 14, 1] and reference therein.

In contrast to the continuous case, here we work only with an ordered sets of
points constructed from the solutions of eq. 12. Still, all the information
about the index can be recovered from them.

We will use the following lemma:

######  Lemma 2.

Suppose that a quadratic form $\mathcal{Q}$ is defined on finite dimensional
vector space $X$ and a subspace $V\subseteq X$ is given. Denote by
$V^{\perp_{\mathcal{Q}}}$ the $\mathcal{Q}-$orthogonal subspace to $V$ i.e.:

$V^{\perp_{\mathcal{Q}}}:=\\{u\in X:\mathcal{Q}(u,v)=0,\forall v\in V\\}.$

Denote by $ind^{-}Q$ the number of negative eigenvalues of $Q$, then:

$\displaystyle ind^{-}\mathcal{Q}$
$\displaystyle=ind^{-}\mathcal{Q}|_{V}+ind^{-}\mathcal{Q}|_{V^{\perp_{\mathcal{Q}}}}$
$\displaystyle+\dim\big{(}(V\cap V^{\perp_{\mathcal{Q}}})/(V\cap
V^{\perp_{\mathcal{Q}}}\cap\ker\mathcal{Q})\big{)}.$

We will apply iteratively this formula to the subspaces $V_{k}:=U_{k}\cap\ker
E^{k+1}$. The subspaces $U_{k}$ were defined in (9) and provide a filtration
of our space of variations. One has a natural filtration $V_{k-1}\subset
V_{k}$ as well. Indeed, if $u\in V_{k-1}$, then $x_{k}=E^{k}(u)=0$. Hence, if
we extend the control $u$ to $V_{k}$ by taking $u_{k}=0$, we get from (8)

$x_{k+1}=E^{k+1}(u)=P^{k+1}_{k}E^{k}(u)=0.$

Notice also that for any $k$, the subspace $V_{k-1}$ has codimension $n$ in
$V_{k}$.

We denote for brevity

$\mathcal{Q}_{k}:=\mathcal{J}^{k}|_{\ker E^{k+1}},\quad
V_{k}^{\perp}:=V_{k}^{\perp_{\mathcal{Q}_{k+1}}},\quad\mathcal{Q}_{k}^{\perp}:=\mathcal{Q}_{k}|_{V_{k-1}^{\perp}}.$

Let us apply Lemma 2 with $\mathcal{Q}=\mathcal{Q}_{k+1}$ and
$V=V_{k}\subseteq V_{k+1}=X$. This gives:

$\displaystyle ind^{-}\mathcal{Q}_{k+1}$
$\displaystyle=ind^{-}\mathcal{Q}_{k}+ind^{-}\mathcal{Q}_{k+1}^{\perp}$
$\displaystyle+\dim\big{(}(V_{k}\cap V_{k}^{\perp})/(V_{k}\cap
V_{k}^{\perp}\cap\ker\mathcal{Q}_{k+1})\big{)}.$

Thus, iteration of this formula gives:

$\displaystyle ind^{-}\mathcal{J}|_{\ker E^{N+1}}$
$\displaystyle=\sum_{k=0}^{N-1}ind^{-}\mathcal{Q}_{k+1}^{\perp}$
$\displaystyle+\dim\big{(}(V_{k}\cap V_{k}^{\perp})/(V_{k}\cap
V_{k}^{\perp}\cap\ker\mathcal{Q}_{k+1})\big{)}.$

What is left to do is to express every term using the discrete Hamiltonian
system given in eq. 12.

First of all, let us describe the subspaces $V_{k}^{\perp}\cap V_{k}$ and
$V_{k}^{\perp}\cap V_{k}\cap\ker\mathcal{Q}_{k+1}$.

######  Lemma 3.

Let $V_{k}$ be defined as above and let $\Pi$ be the vertical subspace

$\Pi=\\{(\lambda,x)\in\mathbb{R}^{2n}\,:\,x=0\\}.$

Then we have the following identifications:

$\begin{split}V_{k}^{\perp}\cap
V_{k}&=\ker(\mathcal{Q}_{k})=\Phi_{k+1}(\Pi)\cap\Pi;\\\ V_{k}^{\perp}\cap
V_{k}\cap\ker(\mathcal{Q}_{k+1})&=\Phi_{k+1}(\Pi)\cap
M_{k+1}^{-1}(\Pi)\cap\Pi.\end{split}$

###### Proof.

A direct consequence of the definitions of $V_{k}$ and $V_{k}^{\perp}$ is that
$V_{k}^{\perp}\cap V_{k}=\ker(\mathcal{Q}_{k})$. Lemma 1 identifies the
elements of the kernel with solutions of the equation (12) with
$x_{0}=x_{k+1}=0$. In terms of the corresponding flow $\Phi_{k+1}$ we get:

$\ker(\mathcal{Q}_{k})=\Phi_{k+1}(\Pi)\cap\Pi.$

Similarly we have $\ker(\mathcal{Q}_{k+1})=\Phi_{k+2}(\Pi)\cap\Pi$. Applying
$M_{k+1}^{-1}$ gives

$\ker(\mathcal{Q}_{k+1})=\Phi_{k+1}(\Pi)\cap M_{k+1}^{-1}(\Pi).$

Combining this with the previous formula proves the second isomorphism in the
statement. ∎

Notice that the previous lemma holds for general LQR. In this particular case,
$M^{-1}_{k+1}(\Pi)$ is transversal to the fibre. Indeed, from the explicit
expression of $M_{k+1}$ we have that $M_{k+1}(\Pi)\cap\Pi\neq 0$ if and only
if $\ker\Gamma_{k+1}(P^{k+1}_{k+2})^{t}\neq 0$, which is never the case under
1. Thus

$\begin{split}V_{k}^{\perp}\cap
V_{k}=\ker(\mathcal{Q}_{k})&=\Phi_{k+1}(\Pi)\cap\Pi;\\\ V_{k}^{\perp}\cap
V_{k}\cap\ker(\mathcal{Q}_{k+1})&=0.\end{split}$

The next step is to compute $ind^{-}\mathcal{Q}_{k}^{\perp}$. To do so
consider the standard symplectic form on $\mathbb{R}^{2n}$:

$\sigma((\lambda_{1},x_{1}),(\lambda_{2},x_{2}))=\lambda_{1}^{t}x_{2}-\lambda_{2}^{t}x_{1}=\begin{pmatrix}\lambda_{1}^{t}&x_{1}^{t}\end{pmatrix}\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}\begin{pmatrix}\lambda_{2}\\\ x_{2}\end{pmatrix}.$

A plane $\Lambda\subset\mathbb{R}^{2n}$ is called Lagrangian if
$\sigma|_{\Lambda}=0$ and $\dim\Lambda=n$, i.e., it is a maximally isotropic
space.

The Maslov form $m(\Lambda_{1},\Lambda_{2},\Lambda_{3})$ of a triple of
Lagrangian spaces $\Lambda_{1},\Lambda_{2},\Lambda_{3}$ is a symmetric
quadratic form on $(\Lambda_{1}+\Lambda_{3})\cap\Lambda_{2}$ defined as
follows. Let $p_{2}=p_{1}+p_{3}$ where $p_{i}\in\Lambda_{i}$. Then

$m(\Lambda_{1},\Lambda_{2},\Lambda_{3})(p_{2})=\sigma(p_{1},p_{3}).$

######  Lemma 4.

The negative index of $\mathcal{Q}_{k}$ restricted to $V_{k-1}^{\perp}$
coincides with the negative index of $m(\Phi_{k}(\Pi),M_{k}^{-1}(\Pi),\Pi)$.

###### Proof.

An element $u$ belongs to $V_{k}$ if and only if it is of the form:

$u=(u_{0},\dots,u_{k-1},F_{k}(u)),$

where $F_{k}$ is defined in (10). Hence we can identify $V_{k}$ with
$\mathbb{R}^{nk}$.

It follows that an element $u$ belongs to $V_{k-1}^{\perp}$ if and only if for
all $v\in V_{k-1}\simeq\mathbb{R}^{n(k-1)}$:

$\displaystyle 0=\mathcal{Q}_{k}(u,v)=\frac{1}{2}\left(\sum_{i=0}^{k}\langle
u_{i},v_{i}\rangle-\sum_{i=1}^{k}\langle Q_{i}E^{i}u,E^{i}v\rangle\right)$
$\displaystyle=\frac{1}{2}\Bigg{(}\sum_{i=0}^{k-2}\langle
u_{i},v_{i}\rangle-\langle(E^{i+1})^{t}Q_{i+1}E^{i+1}u,v\rangle_{L^{2}}+\langle
u_{k-1},F_{k-1}(v)\rangle\Bigg{)}$ $\displaystyle=\frac{1}{2}\langle
u-\sum_{i=0}^{k-2}(E^{i+1})^{t}Q_{i+1}E^{i+1}u+F_{k-1}^{t}(u_{k-1}),v\rangle_{L^{2}}$

where the $L^{2}$-scalar product is defined on functions on a finite set with
values in $\mathbb{R}^{n}$. Here we used that $E^{k}(v)=0$ since $v\in
V_{k-1}$.

We have thus shown that if $u\in V_{k-1}^{\perp}$ then

$\langle
u-\sum_{i=0}^{k-2}(E^{i+1})^{t}Q_{i+1}E^{i+1}u,v\rangle_{L^{2}}=-\langle
F_{k-1}^{t}(u_{k-1}),v\rangle_{L^{2}},$

for all functions $v:\\{0,\dots,k-1\\}\to\mathbb{R}^{n}$.

Let us now simplify the form $\mathcal{Q}_{k}^{\perp}$ using this expression.
Take an element $u\in V_{k-1}^{\perp}$. Then

$\begin{split}&\mathcal{Q}_{k}(u)=\frac{1}{2}\Bigg{(}\left(\sum_{i=0}^{k-2}\langle
u_{i},u_{i}\rangle-\langle u,(E^{i+1})^{t}Q_{i+1}E^{i+1}(u)\rangle\right)\\\
&+\langle u_{k-1},u_{k-1}\rangle+\langle F_{k}(u),F_{k}(u)\rangle-\langle
E^{k}(u),Q_{k}E^{k}(u)\rangle\Bigg{)}\\\ &=\frac{1}{2}(-\langle
F_{k-1}(u),u_{k-1}\rangle+\langle u_{k-1},u_{k-1}\rangle+\langle
F_{k}(u),F_{k}(u)\rangle\\\ &-\langle x_{k},Q_{k}x_{k}\rangle).\end{split}$

Where we used the characterization of $V_{k-1}^{\perp}$ given in the previous
lemma and substituted $E^{k}(u)$ with $x_{k}$. We now rewrite each term in
symplectic form. In particular recall the following relations, which are
consequences of (17) and (3)

$\displaystyle u_{k}$ $\displaystyle=B_{k}^{t}\lambda_{k+1},$ $\displaystyle
x_{k}$ $\displaystyle=P_{k-1}^{k}x_{k-1}+\Gamma_{k-1}\lambda_{k},$
$\displaystyle B_{k-1}^{t}\lambda_{k}$
$\displaystyle=F_{k-1}(u)=-B_{k-1}^{-1}P_{k-1}^{k}x_{k-1}.$

This implies that:

$\displaystyle\langle u_{k-1},u_{k-1}\rangle$
$\displaystyle=\langle\lambda_{k},\Gamma_{k-1}\lambda_{k}\rangle,$
$\displaystyle\langle F_{k-1}(u),u_{k-1}\rangle$
$\displaystyle=-\langle\lambda_{k},P_{k-1}^{k}x_{k-1}\rangle.$

Similarly, $F_{k}(u)$ and $x_{k}$ can be expressed in terms of
$\lambda_{k+1}$.

$-P_{k}^{k+1}x_{k}=B_{k}F_{k}(u)=\Gamma_{k}\lambda_{k+1}.$

We thus get the following expression for $u\in V_{k-1}^{\perp}$ after
substituting the previous relations and grouping separately terms with
$\lambda_{k}$ and $\lambda_{k+1}$:

$\begin{split}\mathcal{Q}_{k}(u)&=\frac{1}{2}(\langle\lambda_{k},P_{k-1}^{k}x_{k-1}+\Gamma_{k-1}\lambda_{k}\rangle\\\
&\quad+\langle\lambda_{k+1},\Gamma_{k}\left(\lambda_{k+1}+(P_{k+1}^{k})^{t}Q_{k}x_{k}\right)\rangle)\\\
&=\frac{1}{2}\left(\langle\lambda_{k},x_{k}\rangle+\langle\lambda_{k+1},\Gamma_{k}\left(\lambda_{k+1}+(P_{k+1}^{k})^{t}Q_{k}x_{k}\right)\rangle\right).\end{split}$

On the other hand let us consider the Maslov form
$m(\Phi_{k}(\Pi),M^{-1}_{k}(\Pi),\Pi)$ and write down its explicit expression.
We have on $M^{-1}_{k}(\Pi)\cap(\Phi_{k}(\Pi)+\Pi)$.

$\displaystyle M_{k}^{-1}\begin{pmatrix}\lambda_{k+1}\\\ 0\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\left((P_{k}^{k+1})^{t}-Q_{k}P_{k+1}^{k}\Gamma_{k}\right)\lambda_{k+1}\\\
-P_{k+1}^{k}\Gamma_{k}\lambda_{k+1}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\lambda_{k}+\nu\\\
x_{k}\end{pmatrix}\in\Pi+\Phi_{k}(\Pi).$ (19)

Thus the value of
$m(\lambda_{k+1})=\sigma((\lambda_{k},x_{k}),(\nu,0))=-\langle\nu,x_{k}\rangle$
is determined by:

$\begin{split}m(\lambda_{k+1})&=-\langle\left((P_{k}^{k+1})^{t}-Q_{k}P_{k+1}^{k}\Gamma_{k}\right)\lambda_{k+1}-\lambda_{k},x_{k}\rangle\\\
&=\langle\lambda_{k+1},\Gamma_{k}(\lambda_{k+1}+(P_{k+1}^{k})^{t}Q_{k}x_{k})\rangle+\langle\lambda_{k},x_{k}\rangle\\\
&=2\mathcal{Q}_{k}(u).\end{split}$ (20)

Here, in the second equality we substituted $\nu$ from (3) and in the third
equality $x_{k}$ from the same equation. It follows that negative inertia
index of $\mathcal{Q}_{k}^{\perp}$ is the same as the negative inertia index
of $m(\Phi_{k}(\Pi),M^{-1}_{k}(\Pi),\Pi)$. ∎

To arrive at formula (6) we need to write down explicitly the kernel of the
Maslov form in terms of the initial covector $\lambda_{0}$. We have

$\begin{pmatrix}\lambda_{k}\\\
x_{k}\end{pmatrix}=\begin{pmatrix}\Phi_{k}^{1}\lambda_{0}\\\
\Phi_{k}^{3}\lambda_{0}\end{pmatrix},\quad\Phi_{k}=\begin{pmatrix}\Phi_{k}^{1}&\Phi_{k}^{2}\\\
\Phi_{k}^{3}&\Phi_{k}^{4}\end{pmatrix}.$

Moreover, $\lambda_{k+1}$ is directly determined by $x_{k}$, see section 3,
and thus by $\lambda_{0}$ as:

$\lambda_{k+1}=-\Gamma_{k}^{-1}P^{k+1}_{k}\Phi_{k}^{3}\lambda_{0}.$

After inserting inside the explicit expression of the Maslov form given in
(20) we find that the following quadratic form is the one to consider:

$\langle\lambda_{0},\Big{(}(\Phi_{k}^{3})^{t}\big{(}(P^{k+1}_{k})^{t}\Gamma_{k}^{-1}(P^{k+1}_{k})-Q_{k}\big{)}\Phi_{k}^{3}+(\Phi_{k}^{3})^{t}\Phi_{k}^{1}\Big{)}\lambda_{0}\rangle.$

Now each term in the formula for $ind^{-}\mathcal{J}|_{\ker E^{N+1}}$ is
identified and collecting all the pieces together finishes the proof of
Theorem 1.

##  4 Proof of the determinant formula

In this section we prove Theorem 2. Consider the quadratic form
$2\mathcal{J}$. As discussed in the introduction, we can pass from control
variables $u$ to state variables $x$. This gives the Jacobi matrix
${\mathcal{I}}$ once we impose the boundary conditions $x_{0}=x_{N+1}=0$. This
matrix decomposes as the sum of two terms. A diagonal matrix $Q$ and a
tridiagonal one $T=(T_{i,j})_{i,j}$ given by:

$\displaystyle{\mathcal{I}}_{s}=T-sQ,\quad
Q=\mathrm{diag}(Q_{i}),\,i=1,\dots,N,$ $\displaystyle
T_{i}=\Gamma_{i-1}^{-1}+(P_{i}^{i+1})^{t}\Gamma^{-1}_{i}P_{i}^{i+1},\quad
T_{i,i+1}=-(P_{i}^{i+1})^{t}\Gamma^{-1}_{i}.$

Notice, moreover, that the matrix $T$ is positive definite, since it is a
multiple of the quadratic form $\sum_{i=0}^{N}|u_{i}|^{2}$. For this reason we
can define its square root $T^{1/2}$ which is again positive definite and
symmetric. We will use this observation to relate the determinant of
$\Phi_{3}(s)$ to the characteristic polynomial of a suitable $nN$ matrix.

######  Lemma 5.

The function $\det(\Phi_{3}(s))$ is a polynomial of degree
$nN-\sum_{k=1}^{N}\dim(\ker(Q_{k}))$ and in particular it is a multiple of
$\det(T-sQ)$.

###### Proof.

The first thing we have to notice is that the entries of $\Phi(s)$ are
polynomials of degree at most $N$ since $\Phi(s)$ is the product of $N$ affine
in $s$ matrices. This implies that $\det(\Phi_{3}(s))$ has at most degree
$nN$.

This is exactly the dimension of the space of controls. We know by Lemma 1
that this polynomial has at least $d:=nN-\sum_{i}\dim(\ker Q_{i})$ roots. This
is because our perturbation $\mathcal{J}_{s}$ sees every eigenvalue of
$\mathcal{J}$ but $1$, whose eigenspace corresponds to the kernel of $Q$. In
particular when the $Q_{i}$ are all invertible, $\det(\Phi_{3}(s))$ has degree
$nN$, vanishes on the same set as $\det(T-sQ)$ and thus they are multiples. In
general we have:

$\begin{split}\det(T-sQ)&=\det(T)\det(1-sT^{-1/2}QT^{-1/2})\\\
&=s^{nN}\det(T)\det(s^{-1}-T^{-1/2}QT^{-1/2})\\\
&=\lambda^{-nN}\det(T)\det(\lambda-T^{-1/2}QT^{-1/2})|_{\lambda=s^{-1}}.\end{split}$

Order the non zero eigenvalues of $T^{-1/2}QT^{-1/2}$ as
$\\{\lambda_{i}\\}_{i=1}^{d}$. Since $\det(\lambda-T^{-1/2}QT^{-1/2})$ has $d$
roots different from zero and a root of order $\sum_{i}\dim(\ker(Q_{i}))$ at
$\lambda=0$ we can rewrite the formula as:

$\det(T-sQ)=\prod_{i=1}^{d}(1-s\lambda_{i})\det(T).$

Thus we have that $\det(\Phi_{3}(s))$ and $\det(T-sQ)$ are both polynomials of
degree $d$ vanishing on the same set and thus they must be multiples. ∎

By Lemma 5 there exists a constant $C\neq 0$ such that

$\det(T-sQ)=C\det\Phi_{3}(s).$

The constant can be computed by evaluating both sides at $s=0$. Using
induction we can prove that

$\Phi(0)=\begin{pmatrix}(P^{0}_{N+1})^{T}&0\\\
\sum_{i=0}^{N}P^{N+1}_{i+1}\Gamma_{i}(P^{0}_{i+1})^{T}&P^{N+1}_{0}\end{pmatrix}.$
(21)

Hence

$\begin{split}\det(\Phi_{3}(s))|_{s=0}&=\det\left(\sum_{i=0}^{N}P^{N+1}_{i+1}\Gamma_{i}(P^{0}_{i+1})^{T}\right)\\\
&=\det\left(\sum_{i=0}^{N}P^{N+1}_{i+1}\Gamma_{i}(P_{i+1}^{N+1})^{T}\right)\det(P^{0}_{N+1}).\end{split}$

Now it only remains to find $\det T$. Recall that $T$ corresponds to the
quadratic form $\frac{1}{2}\sum_{k=0}^{N}|u_{k}|^{2}$.

Given a trajectory of (3) with $x_{0}=x_{N+1}=0$ we can always reconstruct the
control. This gives a bijection between $(u_{0},\dots,u_{N-1})$ and
$(x_{1},\dots x_{N})$ which can be written using a block-triangular matrix

$L=\begin{pmatrix}B_{0}^{-1}&0&\dots&0&0\\\
-B_{1}^{-1}P^{2}_{1}&B_{1}^{-1}&\dots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\dots&-B_{N-1}^{-1}P^{N}_{N-1}&B_{N-1}^{-1}\end{pmatrix}$

Now we wish to use $F_{N}$ (see eq. 10) to pull back to $U_{N-1}$ the
quadratic form:

$\sum_{i=0}^{N}|u_{i}|^{2}$

seen as a quadratic form on $\ker E^{N+1}$. To do so we simply need to resolve
for $u_{N}$ using (11) and plug it inside the form above. This gives the form
$q(v)=\langle(1+F_{N}^{*}F_{N})v,v\rangle$ on $U_{N-1}$.

We can now change coordinates and work with the $x$ variable. One can see that
$T$ can be written as a symmetric matrix on
$\mathbb{R}^{nN}\times\mathbb{R}^{nN}$ of the form:

$T=L^{T}(1+F_{N}^{*}F_{N})L.$

So if we want to compute its determinant, we need to find determinant of $L$
and of $1+F_{N}^{*}F_{N}$. We clearly have

$\det L=\prod_{i=0}^{N-1}B_{i}^{-1}.$

So it only remains to find $\det(1+F_{N}^{*}F_{N})$.

Notice that, since $F_{N}$ is a surjection, it has a large kernel of dimension
$N(n-1)$. If $v\in\ker F_{N}$, then $v$ is an eigenvector of $q$ with
eigenvalue $1$. So the determinant can be computed multiplying the remaining
$n$ eigenvalues. We claim that they are all equal to the eigenvalues of
$1+F_{N}F_{N}^{*}$. Indeed, if $v\notin\ker F_{N}$ is an eigenvalue of $q$,
then

$(1+F_{N}^{*}F_{N})v=\lambda v.$

Applying $F$ to both sides gives

$(1+F_{N}F_{N}^{*})F_{N}v=\lambda F_{N}v.$

Hence we obtain that:

$\begin{split}\det T&=\prod_{i=0}^{N-1}\det(\Gamma_{i})^{-1}\\\
&\times\det\left(1+\sum_{j=0}^{N-1}B_{N}^{-1}P^{N+1}_{j+1}\Gamma_{j}^{-1}(P^{N+1}_{j+1})^{t}(B_{N}^{-1})^{t}\right)\\\
&=\prod_{i=0}^{N}\det(\Gamma_{i})^{-1}\det\left(\Gamma_{N}+\sum_{j=0}^{N-1}P^{N+1}_{j+1}\Gamma_{j}^{-1}(P^{N+1}_{j+1})^{t}\right).\end{split}$

This implies that $C=\det(P_{0}^{N+1})\prod_{k=0}^{N}\Gamma_{k}^{-1}$ and
thus:

$\det(\Phi_{3}(s))\det(P_{0}^{N+1})\prod_{k=0}^{N}\det\Gamma_{k}^{-1}=\det(T-sQ),$
(22)

which gives formula (7) when $s=1$.

##  5 Some examples and applications

###  5.1 Tridiagonal Toeplitz Matrices

Let us compute the determinant of a symmetric tridiagonal Toeplitz matrix

$T_{N}(s,r)=\begin{pmatrix}s&r&0&\dots&0&0\\\ r&s&r&\dots&0&0\\\
0&r&s&\dots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\dots&s&r\\\ 0&0&0&\dots&r&s\end{pmatrix},$ (23)

where $s,r\in\mathbb{R}$ and $r\neq 0$.

Recall that eigenvalues of $T_{N}(s,r)$ are explicitly known [13] and are of
the form

$\lambda_{j}=s+2|r|\cos\frac{\pi j}{N+1},\qquad 1\leq j\leq N.$

So we already a have a way of computing the determinant of the matrix as the
product of all eigenvalues:

$\det T_{N}(s,r)=\prod_{j=1}^{N}\left(s+2|r|\cos\frac{\pi j}{N+1}\right).$
(24)

Let us compute the same determinant using the control theoretic approach.

Using formulas (5) we reformulate the problem of computing $\det T_{N}(r,s)$
as a LQR problem. We choose for $1\leq i\leq N$:

$B_{0}=B_{i}=1,\quad A_{0}=0,\quad A_{i}=-(1+r),\quad Q_{i}=1+r^{2}-s.$

In particular, this immediately implies that

$\Gamma_{i}=1,\qquad P^{N+1}_{0}=(-r)^{N}.$

After comparing with (7), we see that only $\Phi_{3}$ needs to be computed. In
this particular case

$\Phi=\begin{pmatrix}-r^{-1}&r^{-1}(1+r^{2}-s)\\\
-r^{-1}&r^{-1}(1-s)\end{pmatrix}^{N}\begin{pmatrix}1&0\\\ 1&1\end{pmatrix}.$

Denote

$M=\begin{pmatrix}-r^{-1}&r^{-1}(1+r^{2}-s)\\\
-r^{-1}&r^{-1}(1-s)\end{pmatrix}.$

The eigenvalues of this matrix are given by

$\mu_{\pm}=\frac{-s\pm\sqrt{s^{2}-4r^{2}}}{2r}.$ (25)

Assume that $\mu_{+}\neq\mu_{-}$. Then we can diagonalize $M$ using the matrix

$S=\begin{pmatrix}1+r\mu_{-}&1+r\mu_{+}\\\ 1&1\end{pmatrix},$

as can be verified by a direct computation. Thus

$\Phi=S\begin{pmatrix}\mu_{+}^{N}&0\\\
0&\mu_{-}^{N}\end{pmatrix}S^{-1}\begin{pmatrix}1&0\\\ 1&1\end{pmatrix}.$

After all of the matrices are multiplied, one can take the element which lies
under the diagonal and use formula (7). This way, after some simplification,
we obtain an expression which can be compactly written in terms of the two
eigenvalues

$\det
T_{N}(s,r)=(-r)^{N}\frac{(\mu_{+}^{1+N}-\mu_{-}^{1+N})}{(\mu_{+}-\mu_{-})}=(-r)^{N}\sum_{j=0}^{N}\mu_{+}^{j}\mu_{-}^{N-j}.$
(26)

If $\mu_{+}=\mu_{-}=\mu$, then using the continuity of the determinant this
formula reduces to

$\det T_{N}(s,r)=N(-r\mu)^{N}.$

###  5.2 Generalized Euler’s identity

The formula for the determinant of a tridiagonal matrix has been known for a
long time. However, we would like to underline that control theoretic approach
can be conceptually advantageous. To illustrate this claim, we now prove the
generalized Euler’s identity (2).

Notice that we get the classical Euler’s identity if $a=0$:

$\frac{\sin b}{b}=\prod_{j=1}^{\infty}\left(1-\frac{b^{2}}{(\pi
j)^{2}}\right).$

The identity (2) appeared in [2]. Here we give a proof which uses only
elementary tools from linear algebra and analysis. The idea is to take a
certain sequence of tridiagonal Toeplitz matrices of increasing rank and to
compute the asymptotic expansion of the determinants in two different ways:
using formula (24) and (26).

We start with the following family of one-dimensional LQR problems

$\dot{x}=ax+u,\quad\frac{1}{2}\int_{0}^{1}\left(u(t)^{2}-\varepsilon(a^{2}+b^{2})x^{2}(t)\right)dt\to\min.$

and $x(0)=x(1)=0$. Here $\varepsilon\in\\{0,1\\}$. Morally we want to compute
the determinant of the quadratic functional restricted to the curves which
satisfy the constraints. The difficulty here is that the functional in
question is a quadratic form on an infinite-dimensional space. We will compute
the determinant using finite-dimensional approximations coming from
discretizations of this continuous LQR problem with time step equal to $1/N$.
We consider discrete LQR problems of the form

$x_{k+1}-x_{k}=\frac{1}{N}(ax_{k}+u_{k}),\,\,\frac{1}{2N}\sum_{k=0}^{N}u_{k}^{2}-\varepsilon(a^{2}+b^{2})x_{k}^{2}\to\min.$

In order to apply our formula we introduce a new control function
$v=\sqrt{N}u$, which gives an equivalent optimal control problem:

$x_{k+1}-x_{k}=\frac{a}{N}x_{k}+\frac{v_{k}}{\sqrt{N}},\quad\frac{1}{2}\sum_{k=0}^{N}v_{k}^{2}-\varepsilon\frac{a^{2}+b^{2}}{N}x^{2}\to\min.$

Hence we get as parameters of the LQR problem

$\displaystyle A_{k}=aN^{-1},\quad B_{k}=N^{-1/2},$ $\displaystyle
Q_{k}=\varepsilon(a^{2}+b^{2})N^{-1},\quad\Gamma_{k}=N^{-1}.$

We can now use formulas (5) to show that we are actually computing the
determinants of a tridiagonal Toeplitz matrices $T_{N}(r,s)$ from the previous
subsection with parameters:

$r=-(a+N),\qquad
s=N+N\left(1+\frac{a}{N}\right)^{2}-\varepsilon\frac{a^{2}+b^{2}}{N}.$

Let us denote those matrices as $\hat{T}_{\varepsilon}(a,b,N)$.

As discussed previously we want to compute $\det\hat{T}_{\varepsilon}(a,b,N)$
in two different ways and then take the limit $N\to\infty$. Foreshadowing a
bit the computation, the sequence $\det\hat{T}_{\varepsilon}(a,b,N)$ does not
converge. So instead we compute the limit

$\lim_{N\to\infty}\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}.$

There are two motivations for using exactly this ratio. First of all we can
notice from formula (26) that if two systems have the same matrices $A_{k}$
and $B_{k}$, then the ratio of the corresponding determinants depends only on
the associated discrete flows $\Phi$. Since we have started with an
approximation to a continuous system, it is natural to expect that the flows
will converge and this way we obtain a finite quantity. Another reason is the
celebrated Gelfand-Yaglom formula in quantum mechanics [10], which computes
the ratio of determinants of a 1D quantum particle with a potential and a 1D
free quantum particle as a ratio of determinants of two finite-rank matrices.

Let us first use formula (26) to compute the ratio. The two eigenvalues of the
flow matrix now depend on $a,b,N,\varepsilon$. To keep formulas clean we omit
the dependence on $a,b,N$ in our notations and indicate explicit dependence on
$\varepsilon$ through a sub-index or super-index. We also assume that $a>0$,
$b>0$. Other sign choices will be a direct consequence of (2) itself.

First note that $r$ is independent of $\varepsilon$. If we write
$\Delta_{\varepsilon}=s_{\varepsilon}^{2}-r^{2}$, then the formula (25) for
eigenvalues reads as

$\mu^{\varepsilon}_{\pm}=\frac{-s_{\varepsilon}\pm\sqrt{\Delta_{\varepsilon}}}{2r}.$

Formula (26) in this notations reads as

$\displaystyle\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}=\frac{\left((\mu_{+}^{1})^{1+N}-(\mu_{-}^{1})^{1+N}\right)}{\left((\mu_{+}^{0})^{1+N}-(\mu_{-}^{0})^{1+N}\right)}\frac{(\mu_{+}^{0}-\mu_{-}^{0})}{(\mu_{+}^{1}-\mu_{-}^{1})}$
$\displaystyle=\frac{\sqrt{\Delta_{0}}}{\sqrt{\Delta_{1}}}\frac{\left((-s_{1}+\sqrt{\Delta_{1}})^{1+N}-(-s_{1}-\sqrt{\Delta_{1}})^{1+N}\right)}{\left((-s_{0}+\sqrt{\Delta_{0}})^{1+N}-(-s_{0}-\sqrt{\Delta_{0}})^{1+N}\right)}.$

Plugging formulas for $r$ and $s_{\varepsilon}$ we find that

$\Delta_{\varepsilon}=4(1-\varepsilon)a^{2}-4\varepsilon
b^{2}+O(N^{-1}),\qquad N\to\infty.$

Hence

$\lim_{N\to\infty}\frac{\sqrt{\Delta_{0}}}{\sqrt{\Delta_{1}}}=\frac{2a}{2bi}.$

Similarly we find that

$s_{\varepsilon}=2N+O(1),\qquad N\to\infty.$

Hence

$\lim_{N\to\infty}\frac{(-s_{1}\pm\sqrt{\Delta_{1}})^{1+N}}{(2N)^{1+N}}=e^{\pm
ib}$

and

$\lim_{N\to\infty}\frac{(-s_{0}\pm\sqrt{\Delta_{0}})^{1+N}}{(2N)^{1+N}}=e^{\pm
a}.$

Thus after collecting everything together, we find

$\lim_{N\to\infty}\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}=\frac{2a}{2bi}\frac{e^{ib}-e^{-ib}}{e^{a}-e^{-a}}=\frac{a\sin
b}{b\sinh a}$

which is exactly the left side of (2).

Now we use formula (24) to compute the same ratio. We get

$\displaystyle\frac{\prod_{j=1}^{N}\left(N+N(1+\frac{a}{N})^{2}+2(a+N)\cos\frac{\pi
j}{N+1}-\frac{a^{2}+b^{2}}{N}\right)}{\prod_{j=1}^{N}\left(N+N(1+\frac{a}{N})^{2}+2(a+N)\cos\frac{\pi
j}{N+1}\right)}$
$\displaystyle=\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}.$

Expanding and simplifying the denominator gives us

$\displaystyle N+N\left(1+\frac{a}{N}\right)^{2}+2(a+N)\cos\frac{\pi j}{N+1}$
$\displaystyle=\frac{1}{N}\left(a^{2}+4(N^{2}+aN)\cos^{2}\frac{\pi
j}{2(N+1)}\right).$

Hence

$\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}=\prod_{j=1}^{N}\left(1-\frac{a^{2}+b^{2}}{a^{2}+4(N^{2}+aN)\cos^{2}\frac{\pi
j}{2(N+1)}}\right).$

Changing the summation index $j\mapsto N+1-j$ gives

$\frac{\det\hat{T}_{1}(a,b,N)}{\det\hat{T}_{0}(a,b,N)}=\prod_{j=1}^{N}\left(1-\frac{a^{2}+b^{2}}{a^{2}+4(N^{2}+aN)\sin^{2}\frac{\pi
j}{2(N+1)}}\right).$

Since $\sin(y)\leq|y|$, we have the following inequality:

$\frac{4(N^{2}+aN)}{(\pi j)^{2}}\sin^{2}\frac{\pi j}{2(N+1)}\leq
1+\frac{(a-2)N-1}{4(N+1)^{2}}$

Thus for any $\epsilon>0$ there exists $N_{0}\in\mathbb{N}$ such that $\forall
N\geq N_{0}$:

$\frac{a^{2}+b^{2}}{a^{2}+4(N^{2}+aN)\sin^{2}\frac{\pi
j}{2(N+1)}}\geq\frac{a^{2}+b^{2}}{a^{2}+(1+\epsilon)(\pi j)^{2}}$ (27)

Define the following partial product:

$\Pi_{k}^{m}(N)=\prod_{j=k}^{m}\left(1-\frac{a^{2}+b^{2}}{a^{2}+4(N^{2}+aN)\sin^{2}\frac{\pi
j}{2(N+1)}}\right).$

By continuity we have that:

$\displaystyle\lim_{N\mapsto\infty}\Pi_{0}^{N}(N)=\lim_{N\mapsto\infty}\Pi_{0}^{N_{0}}(N)\lim_{N\to\infty}\Pi_{N_{0}}^{N}(N)$
$\displaystyle=\prod_{j=0}^{N_{0}}\left(1-\frac{a^{2}+b^{2}}{a^{2}+(\pi
j)^{2}}\right)\lim_{N\to\infty}\Pi_{N_{0}}^{N}(N).$

Since the second factor is the tail of a convergent product, thanks to the
estimates in eq. 27, it converges to $1$ when we take the limit for
$N_{0}\to\infty$ and get:

$\lim_{N_{0}\to\infty}\prod_{j=0}^{N_{0}}\left(1-\frac{a^{2}+b^{2}}{a^{2}+(\pi
j)^{2}}\right)=\lim_{N\mapsto\infty}\Pi_{0}^{N}(N).$

Thus we have proven (2).

## References

  * [1] A. A. Agrachëv.  Quadratic mappings in geometric control theory.  In Problems in geometry, Vol. 20 (Russian), Itogi Nauki i Tekhniki, pages 111–205. Akad. Nauk SSSR, Vsesoyuz. Inst. Nauchn. i Tekhn. Inform., Moscow, 1988.  Translated in J. Soviet Math. 51 (1990), no. 6, 2667–2734.
  * [2] A. A. Agrachev.  Spectrum of the second variation.  Tr. Mat. Inst. Steklova, 304(Optimal noe Upravlenie i Differentsial nye Uravneniya):32–48, 2019.
  * [3] Andrei Agrachev and Ivan Beschastnyi.  Jacobi fields in optimal control: Morse and Maslov indices.  Nonlinear Anal., 214:Paper No. 112608, 47, 2022.
  * [4] Baranzini S. Agrachev A. and I. Beschastnyi.  Index theorems for graph-parametrized optimal control problems.  Nonlinearity, Submitted.
  * [5] Kerstin Ammann and Gerald Teschl.  Relative oscillation theory for jacobi matrices.  arXiv:0810.5648, 2008.
  * [6] Stefano Baranzini.  Functional determinants for the second variation.  submitted to Journal of Differential Equations, 2022.
  * [7] Richard L. Burden, J. Douglas Faires, and Albert C. Reynolds.  Numerical analysis.  Prindle, Weber & Schmidt, Boston, Mass., 1978.
  * [8] Stephan D. Burton.  The determinant and volume of 2-bridge links and alternating 3-braids.  New York J. Math., 24:293–316, 2018.
  * [9] J. J. Duistermaat.  On the Morse index in variational calculus.  Advances in Math., 21(2):173–195, 1976.
  * [10] Gerald Dunne.  Functional determinants in quantum field theory.  Journal of Physics A: Mathematical and Theoretical, 41, 12 2007\.
  * [11] Sondre Tesdal Galtung and Katrin Grunert.  A numerical study of variational discretizations of the Camassa-Holm equation.  BIT, 61(4):1271–1309, 2021.
  * [12] Alfred Hucht.  The square lattice Ising model on the rectangle III: Hankel and Toeplitz determinants.  J. Phys. A, 54(37):Paper No. 375201, 29, 2021.
  * [13] Said Kouachi.  Eigenvalues and eigenvectors of tridiagonal matrices.  Electron. J. Linear Algebra, 15:115–133, 2006.
  * [14] Yiming Long.  Index Theory for Symplectic Paths with Applications.  01 2002.
  * [15] Luca Guido Molinari.  Determinants of block tridiagonal matrices.  Linear Algebra Appl., 429(8-9):2221–2226, 2008.
  * [16] Hermann Schulz-Baldes.  Sturm intersection theory for periodic jacobi matrices and linear hamiltonian systems.  Linear Algebra and its Applications, 436(3):498–515, 2012.
  * [17] Hermann Schulz-Baldes and Liam Urban.  Space versus energy oscillations of prüfer phases for matrix sturm-liouville and jacobi operators.  Electronic Journal of Differential Equations, 2020(76):1–23, 2020\.