Morse index and determinant of block Jacobi matrices via optimal control

Stefano Baranzini ¹¹1Università degli studi di Torino, Via Carlo Alberto 10, Torino, Italy Ivan Beschastnyi ²²2Univesidade 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.