Finite-time Cholesky Factorization

Introduction

The Cholesky factorization is one of the most studied transformations in numerical linear algebra. For a symmetric positive semi-definite matrix $A \in\mathbb{R}^{n \times n}$, it takes the form

where $L \in\mathbb{R}^{n \times n}$ is lower triangular with $L_{ii} \ge0$. This simple transformation has ubiquitous applications in different disciplines of the natural sciences.

The first documented use of this factorization is due to Benoit [3], who published the work of (Commandant) Cholesky in solving linear least-squares problems arising in geodetic networks. Originally developed to solve normal systems of equations, the method later morphed into a solver for positive definite systems

for some $x, b \in\mathbb{R}^{n}$, through the works of Fox, Wilkinson [9], and Turing [16]. The factorization also provides a way to simplify other problems. As pointed out in [9], it can be used to transform generalized eigenvalue problems of the form

where $A$ is symmetric, $B \in\mathbb{R}^{n \times n}$ is positive definite, and $\lambda\in\mathbb{R}$. Let $B = LL^{\mathsf{T}}$ be the Cholesky factorization of $B$, and use the change of variable $x := L^{-\mathsf{T}}y$. We then have

In other words, we manage to turn a generalized eigenvalue problem into an eigenvalue problem corresponding to a symmetric matrix. While this is rarely used today as a direct method, sparse approximations of $L$ are used to precondition iterative methods like LOBPCG [15]. Lastly, the Cholesky factorization readily allows us to calculate the determinant of $A = LL^{\mathsf{T}}$ through

The ubiquity of the Cholesky transformation rests on a combination of computational efficiency, numerical reliability, and efficient implementations. For positive definite matrices, it provides a unique factorization at substantially lower computational cost than an orthogonalization-based method such as Householder QR. For a matrix of size $n \times n$, the Cholesky decomposition requires approximately $n^3/3$ floating-point operations, compared with approximately $4n^3/3$ for Householder QR. For $n = 10^4$, these counts are approximately $3.3 \times10^{11}$ and $1.3 \times10^{12}$, respectively. At a sustained computing capacity of 2 TFlops, these operation counts correspond to approximately $\frac{1}{6}$ seconds for Cholesky and four times longer for QR. Actual runtimes depend on the implementation and hardware, but the difference grows with matrix size and can be significant when the factorization is repeated within an outer-loop application. The appeal of the Cholesky transformation is not efficiency alone. When $A \succ0$, its factorization is unique [2], and its floating-point computation is backward stable [1]. The semidefinite case ($A \succeq0$) can admit infinitely many such factorizations. However, for a permutation matrix $P$, the pivoted matrix $P^{\mathsf{T}}AP$ has a unique Cholesky factorization [2], and the pivoted Cholesky transformation can also be shown to be backward stable [12]. These properties are supported by reliable LAPACK [4] implementations, which use blocked algorithms and BLAS-3 routines to compute the transformation efficiently; a brief summary is given in Table 1. Together, these three factors explain why the Cholesky method is so widely used.

Routine	Description
xPOTRF	Dense positive-definite matrix; blocked algorithm
xPOTF2	Dense positive-definite matrix; unblocked algorithm
xPSTRF	Positive-semidefinite matrix; complete pivoting
xPSTF2	Positive-semidefinite matrix; unblocked complete pivoting

Previous work

Previous work. While pivoted Cholesky is an established method for factorizing positive semi-definite matrices; its implementation is more involved than the canonical Cholesky method. This can be attributed to the additional $O(n^2)$ cost incurred for finding the next pivot; and the memory overhead for reordering the arrays during pivoting. To avoid this, it is common practice to modify the positive-semidefinite (PSD) matrix such that the canonical Cholesky algorithm can be employed while incurring some error due to the modification.

One common modification is diagonal loading. It entails shifting the spectrum of a PSD matrix $A$ by a factor $\mu\in\mathbb{R}_+$ such that $A + \mu I$ is positive definite. This inexpensive modification is popular in spatial statistics (a.k.a adding a nugget), in non-convex optimization (a.k.a Levenberg-Marquardt stabilization), in inverse problems (a.k.a Tikhonov regularization) and in machine learning (a.k.a ridge regression). It is evident that the estimated Cholesky factor depends non-linearly on $\mu$. Therefore, finding an optimal shift parameter usually requires the specification of additional optimality criterion; whose choice is not unique.

Another family of modifications is referred to as modified Cholesky methods. They originated in non-convex optimization where Newton and quasi-Newton methods required a positive definite approximation of the Hessian at a given iteration. Gill, Murray and Wright proposed strategies that entails finding a diagonal matrix $E \in\mathbb{R}^{n \times n}$ such that $A + E$ is positive definite, well-conditioned and sufficiently close to $A$. Schnabel and Eskow proposed another method that uses Gershgorin bounds. Cheng and Higham developed a different modification using a block $LDL^{\mathsf{T}}$ factorization. A systematic catalogue of modified Cholesky methods can be found in []. Unlike diagonal loading, modified Cholesky methods find little use outside of optimization and numerical linear algebra.

Orthogonal to these modification approaches, there exists another class of methods that seek to eliminate the null space from a PSD matrix. These so-called low-rank methods approximate a matrix as $A \approx L_r L_r^{\mathsf{T}}$, where $r \ll n$ is the rank of the approximation. The simplest way to do this is to truncate the eigendecomposition of $A$. This is more expensive than the Cholesky factorization itself, so it is seldom used. On the other hand, Nyström-based approaches, such as the randomly pivoted Cholesky method, scale as $\mathcal{O}(r^2 n)$ and work extremely well when $A$ has a low-rank structure. Furthermore, this family of randomized methods is also stable for downstream applications involving approximations of different matrix functions of these PSD matrices.

Alternatively, one can mimic a pivoted Cholesky method through deflation. Recall that a canonical Cholesky algorithm breaks down when the pivot, say $\alpha\leq0$. In deflation, we ignore the update on this iteration; setting the corresponding column to zero. This method was primarily used to obtain an approximate Cholesky factor for the purpose of preconditioning linear systems arising from discretizations of PDEs. As noted in [], the decision of when to deflate and when not to remains less understood, therefore deflation remains an approximate Cholesky factorization method.

Lastly, Edelman and Holy recently introduced a new method for factoring semidefinite and indefinite matrices using Jordan algebras. This method has the same structure of a canonical Cholesky algorithm except they use the hemiplex number system. A hemiplex number $z$ is of the form $z = a + m\mu+ n\nu$ where $a, m, n \in\mathbb{R}$ and $1, \mu, \nu$ are the basis elements. The number system has additional axioms as specified in []. In computing the Cholesky factorization of the following matrix:

one requires among other conditions that $a^2 = 0$ and $ab = 1$. This is impossible when $a, b \in\mathbb{R}$. On the other hand, the axioms of hemiplex numbers allow for this. Their implementations suggest that their method can be competitive with pivoted Cholesky and $LDL^{\mathsf{T}}$ decompositions. However, a comprehensive benchmark on practical examples remains wanting.

Taken together, these approaches provide different ways to handle positive semidefinite matrices, either by modifying the matrix, constructing a low-rank approximation, handling vanishing pivots, or extending the underlying number system. This review has focused on approaches relevant to the factorization of dense PSD matrices, rather than providing an exhaustive account of sparse and scalable methods or their applications. For a general introduction to Cholesky factorization, its computation and numerical stability, we refer the reader to Higham (2009).

Contributions

This work introduces a new method for calculating the Cholesky factorization of positive-semi definite matrices without pivoting. Its main contributions are as follows:

1. We formulate Cholesky factorization as a variational problem; and derive dynamical systems whose stationary states provide the required optimality conditions. We derive analytical solutions to these dynamical systems and delineate how these solutions can be used to modify the canonical Cholesky algorithm and its blocked counterpart. We refer to the proposed modifications as finite-time Cholesky methods.

2. We analyze the accuracy of the finite-time Cholesky map and discuss its dependence on the evolution time (of the dynamical system) $T$.

3. We demonstrate with synthetic and practical examples that our methods is competitive with (the blocked and unblocked versions of) the pivoted Cholesky methods.

The remainder of the paper is organized as follows. The finite-time Cholesky algorithm is derived in subsection 2.1 and subsection 2.2. Its convergence rates with respect to the evolution time $T$ is studied in subsection 2.3. Using synthetic matrices with different spectral shapes, the dependence of the accuracy of the method on $T$ is demonstrated in section 3. Lastly, we compare our method with pivoted Cholesky methods on different practical matrices in subsection 3.1.

Finite time Cholesky factorization

Let the space of positive-semi definite matrices be denoted by:

Similarly, the space of lower triangular matrices is denoted by:

Recall that the latter is a linear space whereas the former is a cone in $\mathbb{R}^{n \times n}$. Any $A \in\mathcal{A}$ has a unique factorization $A=LL^{T}$ where $L \in\mathcal{L}$. This is referred to as the Cholesky factorization of $A$. For convenience, we will denote this using a map:

where $\operatorname{chol}:\mathcal{A}\mapsto\mathcal{L}$ is a bijective map for all $A \succ0$. On the other hand, when $A \succeq0$, the factorization is not unique and one usually considers an equivalence class of Cholesky factors:

Variational formulation

The Cholesky factorization of a given $A \in\mathcal{A}$ can be interpreted as a structured variational problem:

Here, $\lVert\cdot\rVert_{F}$ is the Frobenius norm of a matrix. Denoting the matrices $A,L$ in block form such as:

where $\alpha:=A_{1,1}$, $r:=A_{2:n,1}$, $l:=L_{1,1}$, $u:=L_{2:n,1}$ and $\widehat{L}:=L_{2:n,2:n}$; one can show that the cost function in (7) can be decoupled. That is:

This dependence of $J$ on $\widehat{L}$ can be further dropped by examining its factorization in terms of its Schur complement.

We will denote one of the Schur complements by $S := B-\frac{1}{\alpha}rr^{T}$. Since, $A \succeq0 \Rightarrow S \succeq0$ and $uu^{T}=\frac{1}{\alpha}rr^{T}$, we can show that $\widehat{L}\widehat{L}^{T}=S$ which should be recursively computed by solving the variational problem $\mathop{\arg\min}_{\widehat{L}\in\mathcal{L}_{n-1}}\|S-\widehat{L}\widehat{L}^{T}\|_{F}^{2}$. Hence,

The minimizer to this reduced functional can be obtained by solving the zeros of:

When $\alpha>0$ and using the constraint on $\mathcal{L}$ that $l\geq0$, the minimizers are obtained by substitution.

Substituting back, one arrives at the familiar recursive formulation of the Cholesky factorization.

For positive definite matrices $A\succ0$, this procedure for estimating $L$ is backwards stable. However, positive-semi-definite matrices require additional pivoting before chol can be used reliably.

Cholesky factorization as a gradient flow

Finite-time Cholesky factorization does not resort to pivoting. Instead, it solves the local minimization problem using its corresponding normalized gradient flow.

Here, $l_{0}\geq0$ is chosen appropriately (as we will specify later). Since $u(0)=0$ and $\lim_{t\to\infty}u(t)=\frac{1}{\sqrt{\alpha}}r$ for $\alpha>0$ and $\lim_{t\to\infty}u(t)=0$ when $\alpha=0$, it is sensible to constrain the dynamics of $u$ along the direction of the residual. Henceforth, we assume that, there exists $g:[0,T]\to\mathbb{R}$ such that $g(0)=0$ and $u(t)=rg(t)$. The corresponding reduced system of ODEs is given by:

To the best of our knowledge, there does not exist an analytical solution to this coupled ODE. As solving ODEs using numerical integrators for each iteration in a Cholesky decomposition function call is expensive, we will instead consider a surrogate ODE that has an analytical solution.

Its solution is given by:

A proof of this derivation is found in the Appendix. Using this, the finite-time Cholesky factorization of the matrix $A$ in (8) is given by:

Notice that this surrogate ODE has the same stationary points as the reduced ODE. The only caveat of deleting the feedback term from its vector field is that $g(t)l(t)=1$ no longer holds. We can argue that this effect is not really necessary. For one, when $A \succ0$, $\lim_{t\to\infty}g(t)l(t)=1$. Furthermore when $A \succeq0$, $r^{T}r=0$. Therefore, for computing the coefficients in the Cholesky decomposition, the surrogate ODE should suffice.

The issue of how the terminal time $t=T$ is chosen has to be taken up on a case-by-case basis. The initial condition $l_0$ on the other hand is simply chosen as $l_0=\max\{1,\alpha^2\}$.

The entire finite-time Cholesky factorization algorithm is specified in Algorithm 2.1. $\blacksquare$

A blocked implementation of the same algorithm is listed in Algorithm 2.2. As the two implementations are equivalent in infinite precision arithmetic, it suffices to analyze the properties of the former.

Algorithm 2.1 FT-Cholesky factorization

Require: Symmetric PSD matrix $A \in\mathbb{R}^{n \times n}$, finite horizon $T > 0$.

Ensure: Lower-triangular finite-time factor $L_T$

1. Define

2. $L \leftarrow0$

3. for $k = 1,\ldots,n$ do

4. $\alpha_k \leftarrow A_{kk} - L_{k,1:k-1}L_{k,1:k-1}^{\mathsf T}$

5. $r_k \leftarrow A_{k+1:n,k} - L_{k+1:n,1:k-1}L_{k,1:k-1}^{\mathsf T}$

6. $s_k \leftarrow\max\{1,\alpha_k\}$

7. $p_k \leftarrow\left[s_k^{-1}e^{-2\alpha_kT} + 2T\varphi_1(-2\alpha_kT)\right]^{-1}$

8. $\ell_k \leftarrow\sqrt{p_k}$

9. $g_k \leftarrow T\ell_k\varphi_1(-p_kT)$

10. $L_{kk} \leftarrow\ell_k$

11. $L_{k+1:n,k} \leftarrow g_kr_k$

12. end for

13. $L_T \leftarrow L$

14. return $L_T$

Algorithm 2.1 has a computational complexity of $O(n^3/3)$. In comparison to the canonical Cholesky method, the method requires evaluation of two additional scalar functions, namely $x \mapsto e^x$ and $x \mapsto(e^x-1)/x$ per iteration. These are standard functions in the GNU scientific library [10] and can be included in existing Cholesky implementations without adding much computational overhead. Algorithm 2.2 has the same computational complexity. However, it reorganizes Algorithm 2.1 to use BLAS-3 kernels [8]. Therefore, it utilizes the underlying hardware more efficiently, particular when $n$ is large.

We will now demonstrate these algorithms on two positive-semi-definite matrix systems, namely the Hilbert matrices (denoted by $H \in\mathbb{R}^{n \times n}$) and the periodic finite-difference 1D-Laplacian matrices $D_{\Delta} \in\mathbb{R}^{n \times n}$. They are defined as follows:

The former is (strictly positive definite but) notoriously ill-conditioned for large $n$, and the latter satisfies $D_{\Delta}\mathbf{1}=0$. Furthermore, in finite precision arithmetic, the representations of these matrices can also have negative eigenvalues. Promptly, one would expect canonical Cholesky routines like xPOTF2 and xPOTRF to terminate prematurely. The finite time routines, on the other hand, return Cholesky factors that accurately reconstruct the matrices. This is shown in Table 2.

Method	Hilbert matrix	Periodic Laplacian
FT-Cholesky, unblocked	$1.69 \times10^{-14}$	$1.10 \times10^{-16}$
FT-Cholesky, blocked	$1.65 \times10^{-14}$	$1.10 \times10^{-16}$
DPOTRF	—	—
DPOTF2	—	—

Algorithm 2.2 Block FT-Cholesky factorization

Require: Symmetric PSD matrix $A \in\mathbb{R}^{n \times n}$, finite horizon $T > 0$, block size $b \ge1$.

Ensure: Lower-triangular finite-time factor $L_T$

1. Define

   $$\varphi_1(z) = \begin{cases} \dfrac{e^z-1}{z}, & z \ne0, \\ 1, & z = 0. \end{cases}$$

2. $L \gets0$

3. $S \gets A$

4. for $k = 1, 1+b, \ldots, 1+b\lfloor(n-1)/b\rfloor$ do

5. $q \gets\min(b,n-k+1)$

6. $I \gets k:k+q-1$

7. $J \gets k+q:n$

8. $B \gets S_{I,I}$

9. $C \gets S_{J,I}$

10. $U \gets0_{q \times q}$

11. $V \gets0_{\lvert J\rvert\times q}$

12. for $j = 1,\ldots,q$ do

13. $\alpha_j \gets B_{jj}-U_{j,1:j-1}U_{j,1:j-1}^{\mathsf{T}}$

14. $r_j \gets B_{j+1:q,j}-U_{j+1:q,1:j-1}U_{j,1:j-1}^{\mathsf{T}}$

15. $s_j \gets C_{:,j}-V_{:,1:j-1}U_{j,1:j-1}^{\mathsf{T}}$

16. $\sigma_j \gets\max\{1,\alpha_j\}$

17. $p_j \gets[\sigma_j^{-1}e^{-2\alpha_jT}+2T\varphi_1(-2\alpha_jT)]^{-1}$

18. $\ell_j \gets\sqrt{p_j}$

19. $g_j \gets T\ell_j\varphi_1(-p_jT)$

20. $U_{jj} \gets\ell_j$

21. $U_{j+1:q,j} \gets g_jr_j$

22. $V_{:,j} \gets g_js_j$

23. end for

24. $L_{I,I} \gets U$

25. $L_{J,I} \gets V$

26. $S_{J,J} \gets S_{J,J}-VV^{\mathsf{T}}$

27. end for

28. $L_T \gets L$

29. return $L_T$

Dependence on $T$

Note that in Table 2, the time evolution time $T$ of the finite-time algorithms was chosen arbitrarily as $T = 5 \times10^{13}$. The dependence of the reconstruction accuracy on $T$ is not immediately obvious. This is studied next. We will first consider positive definite matrices.

Proof. We will first establish arguments for the well-defined-ness of the Cholesky map. This is followed by a slight digression that bounds the error induced by the finite-time evolution. Lastly, we will establish the exponential rate of convergence of the finite time Cholesky map to the canonical Cholesky map using induction. These arguments are detailed below.

Step 1: Well-definedness. Recall that the finite time Cholesky map is given by:

To prove well-defined-ness, it suffices to show that the Schur complement: $B-g(T)^2rr^{\mathsf{T}} \succ0$ for all $T > 0$. To prove this, consider the difference between the finite-time Schur complement and its canonical counterpart.

Since, $S \succ0$, it follows that if $\frac{1}{\alpha} > g(T)^2$, then $S_T \succ0$. To show this, regroup the activation functions as:

Since, $l_0^2 \geq\alpha$ by the construction of the FT-Cholesky algorithm, it follows that $1/l(T)^2 \leq1/\alpha$. Recalling that for $\alpha> 0$:

we can conclude that $g(T)^2 < \frac{1}{\alpha}$. Therefore, $L_T$ is well defined for all $T > 0$.

Step 2: Bounds on scalar errors. Bounds for the absolute difference between the scalar finite-time and canonical Cholesky coefficients are derived below. Note that:

Since, $l_0^2 \geq\alpha$; $\phi(\alpha,l_0):=\left|\frac{\alpha}{l_0^2}-1\right|=\left(1-\frac{\alpha}{l_0^2}\right)$. It follows that:

Therefore, $\left|\sqrt{\alpha}-l(T)\right|=\mathcal{O}(e^{-2\alpha T})$. Similarly,

If $l_0^2\geq\alpha$, then $l(T)\geq l_0\geq\sqrt{\alpha}$; it follows that:

As a result, $\left|\frac{1}{\sqrt{\alpha}}-g(T)\right|=\mathcal{O}(e^{-\alpha T})$.

Step 3: Convergence rate. As specified in the proof sketch, we will derive convergence rates using induction. We will deduce the hypothesis for the convergence rates by examining two base cases, namely $n=2,3$. For a matrix given matrix $A_2\succ0\in\mathbb{R}^{2\times2}$, we denote its canonical and finite time Cholesky factors by $[L]_2$ and $[L_T]_2$ respectively. They are given by:

Denoting the Schur complements by $s,s_T$, it follows that:

We may bound the last term as follows:

Using $g_T(a_{11})\leq\frac{1}{\sqrt{a_{11}}}$, $s_T\geq s$ and the previous scalar error bound, it follows that:

Using this, we can show that

Similarly, let $A_3\in\mathbb{R}^{3\times3}$ be symmetric positive definite, and let $[L]_3$ and $[L_T]_3$ denote its canonical and finite-time Cholesky factors, respectively. Define the exact and finite-time Schur-complement entries by

Then, the corresponding Cholesky factors are given by:

Proceeding in a similar fashion as before, we can show that:

Note that: $l_{11}^{2}=a_{11}$, $l_{22}^{2}=a_{22}-\frac{a_{12}^{2}}{a_{11}}$. This constitutes the base case for the induction argument. Assume that for analogously defined positive definite matrix $A_{k-1}$ and its Cholesky factors $[L]_{k-1}$ and $[L_T]_{k-1}$ for some $k-1<n$, the following is true:

Denoting $A_k\in\mathbb{R}^{k\times k}$ and its finite time Cholesky factor $[L_T]_k$ by:

the error in the factorization is given by:

Since the $\operatorname{Chol}:A\mapsto L$ is Lipschitz continuous say with constant $K>0$, it follows that:

We will first bound the first term to the right using the scalar bounds derived for $A_2$.

From step 1, we know that $s_{2:k,2:k,T}\succ0$, $\forall T>0$. Using the induction hypothesis, it follows that:

Here, $\nu_{k-1}$ is the minimum of the square of the diagonal entries of $\operatorname{Chol}(s_{2:k,2:k,T})$. This quantity can in-turn be related to $l_{ii}$ where $i \in\{2,\ldots,k\}$. Firstly, $s_{2:k,2:k,T} \succeq s_{2:k,2:k}$. Therefore, $\operatorname{diag}(\operatorname{Chol}(s_{2:k,2:k,T}))_j \ge\operatorname{diag}(\operatorname{Chol}(s_{2:k,2:k}))_j$ for $j \in\{1,\ldots,k-1\}$. Additionally, $\operatorname{diag}(\operatorname{Chol}(s_{2:k,2:k}))_j = l_{j+1,j+1}$. Consequently,

Therefore, for some constants $C_1,C_2>0$:

This completes the inductive argument. Hence the result is proven. $\square$

This convergence rate is verified numerically for a Matern-5/2 matrix in Figure 1. We will now turn our attention to positive semi-definite matrices. Since, Cholesky factorization for such matrices over $\mathbb{R}$ is not unique, we will settle for the modest goal of bounding the growth of the reconstruction error with respect to $T$.

Proof. The proof borrows several elements from and has the same structure as Theorem 2.1. Its arguments are detailed below.

Step 1: Well-definedness. Recall that the finite-time Cholesky map is well-defined if:

Since, $A\succeq0 \Longrightarrow S\succeq0$, it reduces to verifying $1/\alpha\ge g(T)^2$. Here, there are two cases. When $\alpha>0$, we know from the proof of Theorem 2.1 that $S_T\succ0$. On the other hand, when $\alpha= 0 \Longrightarrow r = 0$ and $g(T) < \infty$. Therefore, $S_T = B$ which is positive semi-definite by definition.

Step 2: Bounds on scalar errors. The activation functions $l, g$ have different asymptotic behavior as $\alpha\to0$. We can show that:

Using this, we can bound the scalar error functions appearing in Step 3. Firstly, $\forall\alpha> 0$:

Since, $l_0^2 \ge\alpha$ by construction; the denominator of this expression $\ge1$, therefore:

We are also interested in the case $\alpha= 0$.

Secondly, it will be of interest to obtain the bounds for $e^{-Tl_T(\alpha)^2}$. Since, $l_T^2(\alpha)\leq\alpha$, it follows that $e^{-l_T^2(\alpha)T}\leq e^{-\alpha T}$. In other words, $e^{-l_T(\alpha)^2T}=\mathcal{O}(e^{-\alpha T})$.

Step 3: Convergence rate. The reconstruction error rate of $A$ will be obtained by induction. We use the definition in (21). Firstly, we observe that for $A_2\in\mathbb{R}^{2\times2}\succeq0$,

Furthermore when $s_T>0$,

Here, we slightly modify the algorithm and maintain a running update of $l_0$. In this case $l_0$ is $\max\{1,s_T,a_{11}\}$. Since $e^{2s_TT}-1>2s_TT$, it follows that:

Similarly, when $s_T$ approaches zero:

Assuming that the modified $l_0$ is used in both evaluations of the activation function $l$, it follows that:

Similarly, for $A_3 \succeq0$ and its corresponding Cholesky factor $[L_T]_3$, we can use the definitions in Proof 1, proceed in an analogous fashion and show that:

Here, $l_0=\max\{1,a_{11},s_{22,T},s_{33,T}-g_T^2(s_{22,T})s_{23,T}^2\}$. The last term is dependent on $T$ and can be bounded using the following:

Therefore, $\lVert A_3-[L_T]_3[L_T]_3^{\mathsf{T}}\rVert_F=O(T^{-1})$. These two examples constitute the base cases for the induction hypothesis. For $A_k\succeq0$, we hypothesize that:

Then, the proposition can be proven for $A_{k+1}\succeq0$ as follows. We first define:

Here, $\lambda_T=(L_T)_{k+1,1:k}$ and $a\in\mathbb{R}^k$. Then,

Denoting by $r_T:=a_{k+1,k+1}-\lambda_T^{\mathsf{T}}\lambda_T$ and choosing $l_0\geq a_{k+1,k+1}$ note that:

Lastly, it suffices to show that $\lVert a-[L_T]_k\lambda_T\rVert^2$ is bounded.

To this end, we use the fact that $r_T\geq0$ to show that

Furthermore, defining by $x_i:=a_{ii}-\sum_{j=1}^{i-1}([L_T]_k)_{ij}^2$ we can show that

Taken together, we have shown that $\lVert a-[L_T]_k\lambda_T\rVert^2=O(T^{-1})$. Using this and the induction hypothesis $\lVert A_{k+1}-[L_T]_{k+1}[L_T]_{k+1}^{\mathsf{T}}\rVert=O(T^{-1})$.

This rate is verified for a Matern-covariance matrix with controlled null-space in Figure 2. Hence, we have shown through Theorem 2.1 and Theorem 2.2 that for positive semi-definite matrices the accuracy of the finite-time Cholesky algorithm increases monotonously in the evolution parameter $T$ so long as it is well-defined.

Numerical experiments

We will now use the finite-time Cholesky algorithm and its blocked implementation to factor a collection of synthetic matrices. These matrices are randomly generated and are designed to have different conditioning and spectral profiles. They are constructed as follows. Let $W \in\mathbb{R}^{n \times n}$ be a random matrix whose elements $W_{ij} \sim\mathcal{N}(0,1)$. The column space corresponding to $W$ is orthogonalized using a QR decomposition, such that $W = QR$. The synthetic matrix is defined as:

Here, $\Lambda_s \in\mathbb{R}^{n \times n}$ is a diagonal matrix. Its smallest and largest non-zero elements say $[\Lambda_s]_{\min}, [\Lambda_s]_{\max}$ specify the condition number $\kappa$. Furthermore, $[\Lambda_s]_{ii} \in[0,\infty)$ where $i \in[1,n]$. The distribution of these elements on $n$ points on this interval, in turn, specifies the spectral shape of $A_s$. The spectral distributions considered in this work are specified in Table 3.

Spectrum $s$	Eigenvalue $\lambda_i^{(s)}$	Index range
Exponential	$a^{1-(i-1)/(n-1)}$	$i=1,\ldots,n$
Polynomial	$(i/n)^p$	$i=1,\ldots,n$
Log-uniform	$a^{1-u_i}$	$i=1,\ldots,n$
Uniform	$a+(1-a)u_i$	$i=1,\ldots,n$
One-cluster	$a$	$i=1$
	$0.9+0.1(i-2)/(n-2)$	$i=2,\ldots,n$
Two-cluster	$a(1+0.1(i-1)/(m-1))$	$i=1,\ldots,m$
	$0.9+0.1(i-m-1)/(m-1)$	$i=m+1,\ldots,n$

Our experiment entails estimating the finite-time Cholesky factorization, say $L_T$, of a given matrix $A$ and measuring the Frobenius norm relative residual:

We will consider one set of matrices from each of the groups in Table 3. Each of the matrices in a group are specified the condition numbers $\kappa\in\{10^2,10^4,10^6,\ldots,10^{16}\}$. The results of this experiments for one realization of these synthetic matrices is furnished in Figure 3.

Subsequently, we will modify the construction of these synthetic matrices to make them semidefinite. This is realized by randomly assigning certain eigenvalues in the spectrum to zero. We refer to the number of zero eigenvalues of a given matrix as its nullity. The reconstruction errors for nullities 1, 2, 5, 10 for the same synthetic matrices are specified in Figure 4, Figure 5, Figure 6 and Figure 7 respectively.

Observe that the symmetric positive definite set of experiments produces reconstruction errors that are close to the numerical floor of Float64 representations (See Figure 3). On the other hand, matrices with nullities have relatively worse reconstruction errors. The issue is particularly dire for clustered spectral shapes. This can be attributed to the large evolution times (here $T = 10^{14}$) used in these experiments. A modified algorithm that finds the Pareto front over $T$ for which the reconstruction errors are optimal can ameliorate some of these reconstruction issues.

Furthermore, some of these matrices are numerically indefinite in Float64 representations. This will consequently spawn negative pivots in the FT-Cholesky algorithm – a regime where the FT-Cholesky map is not well-posed. Nonetheless, when the negative pivots are in a sufficiently small $\epsilon$ ball around zero, the FT-Cholesky algorithm works without introducing spurious compounding errors.

Applications

In the rest of this section, we will test the FT-Cholesky algorithms on matrices that arise in different applications. Particular attention is given to a litany of kernel matrices and correspondingly to finite-difference matrices of linear elliptic differential operators. We are again interested in the relative Frobenius residual of the factorization. Comparisons are made with xPOTRF, xPSTRF and its blocked implementations.

Kernel matrices

Kernel estimation, RBF interpolation, and Gaussian process regression are related and, under some assumptions, equivalent methods for estimating a function from noisy observations. Let $\mathcal{X} \subset\mathbb{R}^{d}$ and $\mathcal{Y} \subset\mathbb{R}$ denote the input and output spaces, respectively. A target function $f:\mathcal{X} \mapsto\mathcal{Y}$ is observed on certain collocation points $\widehat{\mathcal{X}} \subset\mathcal{X}$, denoted by $\widehat{\mathcal{Y}} := f(\widehat{\mathcal{X}}) + \mathcal{E}$ where $\mathcal{E}$ denotes the appropriate pointwise i.i.d. Gaussian noise with variance $\lambda$. The kernel estimation problem tries to approximate $f$ using an estimator:

Here, $k:\mathcal{X} \times\mathcal{X} \mapsto\mathbb{R}$ is a positive-semi-definite function known as the kernel function. Under certain assumptions, these kernel functions are often Green’s functions of certain linear elliptic differential operators. Common examples of such kernel functions include the Matern kernel, the Wendland kernel, the Gaussian kernel and the like. For a comprehensive review on kernel functions, see [Wendland]. To estimate $g$, one solves the following linear least squares problem:

Furthermore, the average incurred in estimating $g$ with finite information is given by its power function (it’s also the posterior covariance in Gaussian process regression). For $x,y \in\mathcal{X}$, it is defined as:

Both these quantities – $a,\mathcal{P}$ can be easily estimated once the Cholesky factor of $k(\widehat{\mathcal{X}},\widehat{\mathcal{X}})$ is estimated. However, in applications where $\lambda$ is relatively small, the canonical Cholesky implementation terminates prematurely. Therefore, we will carry out this factorization using our finite-time algorithms.

The reconstruction errors of kernel matrices due to Matern, integrated-Wiener process and Wendland priors using the finite-time Cholesky algorithm and its blocked version are furnished in Table 4 and Table 5, respectively. These matrices are constructed using randomly sampled points on an interval of $\mathbb{R}$. The nullities of these matrices are controlled by overlapping points in the collocation set. While this appears contrived, uniform random sampling of collocation points has similar tendencies.

Kernel	Nullity	Finite-time error	DPOTF2	Pivoted rank	Pivoted error
Matérn-5/2	1	$3.629 \times10^{-16}$	Pass	499	$3.241 \times10^{-16}$
	2	$3.596 \times10^{-16}$	Fail	498	$3.253 \times10^{-16}$
	5	$3.643 \times10^{-16}$	Fail	495	$3.216 \times10^{-16}$
	10	$3.664 \times10^{-16}$	Fail	490	$3.253 \times10^{-16}$
IWP-1	1	$5.434 \times10^{-16}$	Pass	499	$3.073 \times10^{-16}$
	2	$5.640 \times10^{-16}$	Fail	498	$3.052 \times10^{-16}$
	5	$5.955 \times10^{-16}$	Fail	495	$3.054 \times10^{-16}$
	10	$6.592 \times10^{-16}$	Fail	490	$3.082 \times10^{-16}$
Wendland	1	$1.764 \times10^{-16}$	Fail	499	$3.515 \times10^{-16}$
	2	$1.824 \times10^{-16}$	Fail	498	$3.500 \times10^{-16}$
	5	$1.949 \times10^{-16}$	Fail	495	$3.514 \times10^{-16}$
	10	$2.229 \times10^{-16}$	Fail	490	$3.552 \times10^{-16}$

Table 5: Reconstruction errors for singular kernel matrices generated by exact repeated points. All matrices have size $n = 500$. Revised blocked FT-Cholesky uses $T = 10^{14}$ and block size $b = 64$. The error is $\lVert K - LL^{\mathsf{T}}\rVert_{\mathrm{F}}/\lVert K\rVert_{\mathrm{F}}$, accounting for permutation for pivoted factors. Standard and pivoted Cholesky use DPOTRF and DPSTRF; pivoted rank uses LAPACK’s default tolerance.

Kernel	Nullity	Finite-time error	DPOTRF	Pivoted rank	Pivoted error
Matérn-5/2	1	$2.284 \times10^{-16}$	Pass	499	$2.113 \times10^{-16}$
	2	$2.301 \times10^{-16}$	Fail	498	$2.140 \times10^{-16}$
	5	$2.343 \times10^{-16}$	Fail	495	$2.129 \times10^{-16}$
	10	$2.400 \times10^{-16}$	Fail	490	$2.142 \times10^{-16}$
IWP-1	1	$5.177 \times10^{-16}$	Pass	499	$1.834 \times10^{-16}$
	2	$5.308 \times10^{-16}$	Fail	498	$1.821 \times10^{-16}$
	5	$5.715 \times10^{-16}$	Fail	495	$1.817 \times10^{-16}$
	10	$6.375 \times10^{-16}$	Fail	490	$1.833 \times10^{-16}$
Wendland	1	$1.567 \times10^{-16}$	Fail	499	$2.791 \times10^{-16}$
	2	$1.657 \times10^{-16}$	Fail	498	$2.825 \times10^{-16}$
	5	$1.819 \times10^{-16}$	Fail	495	$2.801 \times10^{-16}$
	10	$2.109 \times10^{-16}$	Fail	490	$2.799 \times10^{-16}$

Note that both of our methods are competitive with pivoted routines while preserving the memory access patterns of the canonical Cholesky method. While pivoted routines have the advantage of being rank-revealing, it is seldom interesting in practice to estimate the rank of these Gram matrices.

Remark 3.1. Approximation bounds (along with empirical investigations) on the estimating $a,\mathcal{P}$ will be pursued in a subsequent work.

Hermite kernel matrices

A natural generalization of kernel estimation involves approximation, where $f:\mathcal{X} \mapsto\mathcal{Y}$ is smooth, i.e. $C^{k}(\mathcal{X},\mathcal{Y})$. Furthermore, we assume that, in addition to function values on $\widehat{\mathcal{X}}$, we also measure certain derivatives of $f$, namely $Df,D^{2}f,\ldots$ at the same points. We denote the measurement sets by

$\hat{y}^{s}_{k=0} := D^{k}f(\hat{X}) + \mathcal{E}_{k}$ where $\mathcal{E}_{k}$ is again pointwise i.i.d Gaussian noise with variances $\lambda_{k} > 0$. The estimator is specified using the following ansatz:

Here, $[D^{k}]_{1}$ denotes the $k$th derivative operator with respect to the first argument of the kernel $k$, and $b_{j}^{k} \in\mathbb{R}$ denotes the coefficients to be estimated. To estimate these parameters, we consider the derivatives of this ansatz:

Taken together, this estimation procedure with derivative information induces a Gram matrix, which we denote by:

Here, each block say $k(\hat{X},\hat{X}) \in\mathbb{R}^{M \times M}$. Furthermore, each block is not positive semi-definite. On the other hand, $\Theta_{\hat{X},s,M} \succeq0$. We will refer to such structures as Hermite-kernel matrices, eponymous with the Hermite-Birkhoff interpolation problem. An immediate application of this is in solving partial differential equations using collocation. For instance, when solving the Poisson equation:

where $p \in L^2(\Omega)$, one assembles the following Hermite-kernel matrix:

Depending on the smoothness of the kernel, these matrices can be highly ill-conditioned. In the following, we will compare the reconstruction accuracies corresponding to different integer Matern kernels. To control the null space of these matrices, we deliberately introduce duplicate points when the collocation points are sampled. The relative reconstruction accuracies of the unblocked and blocked finite-time Cholesky implementations are furnished in Table 6 and Table 9, respectively.

Table 6: Reconstruction errors for mixed first-order Poisson Hermite Gram matrices with exact repeated operator observations. All matrices have size $n = 1000$. Revised unblocked FT-Cholesky uses $T = 10^{14}$. The error is $\lVert K - LL^\top\rVert_F/\lVert K\rVert_F$, accounting for permutation for pivoted factors. Standard and pivoted Cholesky use DPOTF2 and DPSTF2; pivoted rank uses LAPACK’s default tolerance.

Kernel	Repeats	Finite-time error	DPOTF2	Pivoted rank	Pivoted error
Matérn-3/2	0	$5.157 \times10^{-16}$	Pass	1000	$5.955 \times10^{-16}$
Matérn-3/2	1	$5.166 \times10^{-16}$	Fail	999	$5.875 \times10^{-16}$
Matérn-3/2	2	$5.195 \times10^{-16}$	Fail	998	$5.954 \times10^{-16}$
Matérn-3/2	5	$5.169 \times10^{-16}$	Fail	995	$5.886 \times10^{-16}$
Matérn-3/2	10	$5.248 \times10^{-16}$	Fail	990	$5.847 \times10^{-16}$
Matérn-5/2	0	$5.499 \times10^{-16}$	Pass	1000	$7.752 \times10^{-16}$
Matérn-5/2	1	$5.520 \times10^{-16}$	Fail	999	$7.582 \times10^{-16}$
Matérn-5/2	2	$5.518 \times10^{-16}$	Fail	998	$7.618 \times10^{-16}$
Matérn-5/2	5	$5.483 \times10^{-16}$	Fail	995	$7.592 \times10^{-16}$
Matérn-5/2	10	$5.480 \times10^{-16}$	Fail	990	$7.539 \times10^{-16}$
Matérn-7/2	0	$5.539 \times10^{-16}$	Pass	1000	$6.322 \times10^{-16}$
Matérn-7/2	1	$5.535 \times10^{-16}$	Fail	999	$6.304 \times10^{-16}$
Matérn-7/2	2	$5.561 \times10^{-16}$	Pass	998	$6.320 \times10^{-16}$
Matérn-7/2	5	$5.529 \times10^{-16}$	Fail	995	$6.281 \times10^{-16}$
Matérn-7/2	10	$5.466 \times10^{-16}$	Fail	990	$6.324 \times10^{-16}$

Operator	Boundary	Power	FT error	DPOTRF	Pivoted rank	Pivoted error
Polyharmonic	Neumann	1	$3.653 \times10^{-16}$	Fail	499	$1.508 \times10^{-15}$
		2	$4.486 \times10^{-16}$	Fail	499	$6.552 \times10^{-16}$
		3	$5.549 \times10^{-16}$	Pass	499	$1.361 \times10^{-15}$
Polyharmonic	Periodic	1	$3.793 \times10^{-16}$	Fail	499	$1.514 \times10^{-15}$
		2	$4.111 \times10^{-16}$	Fail	499	$7.078 \times10^{-16}$
		3	$4.150 \times10^{-16}$	Pass	499	$1.510 \times10^{-15}$
Matérn	Neumann	1	$4.056 \times10^{-15}$	Fail	499	$1.802 \times10^{-15}$
		2	$4.000 \times10^{-15}$	Fail	499	$2.843 \times10^{-15}$
		3	$5.662 \times10^{-15}$	Fail	499	$6.049 \times10^{-15}$
Matérn	Periodic	1	$3.978 \times10^{-15}$	Fail	499	$1.805 \times10^{-15}$
		2	$4.183 \times10^{-15}$	Fail	499	$2.869 \times10^{-15}$
		3	$5.886 \times10^{-15}$	Fail	499	$6.069 \times10^{-15}$

We observe that both the unblocked and blocked implementations successfully complete factorization for all the considered kernel matrices. Furthermore, it is competitive with LAPACK’s pivoted Cholesky implementations. The un-pivoted implementations on the other hand, terminate early upon encountering a zero pivot.

Elliptic difference operators

Lastly, we consider finite-difference matrices arising from discretizing second order linear-elliptic differential operators. In free space, these operators are closely related to the positive-definite kernels considered in the previous sections. We consider, for $m \in\mathbb{N}$ and $\nu+ d/2 \in\mathbb{N}$ and $\kappa\in\mathbb{R}$:

Table 7: Reconstruction errors for mixed first-order Poisson Hermite Gram matrices with exact repeated operator observations. All matrices have size $n = 1000$. Revised blocked FT-Cholesky uses $T = 10^{14}$ and block size $b = 64$. The error is $\lVert K - LL^\top\rVert_F/\lVert K\rVert_F$, accounting for permutation for pivoted factors. Standard and pivoted Cholesky use DPOTRF and DPSTRF; pivoted rank uses LAPACK’s default tolerance.

Kernel	Repeats	Finite-time error	DPOTRF	Pivoted rank	Pivoted error
Matérn-3/2	0	$2.068 \times10^{-16}$	Pass	1000	$5.113 \times10^{-16}$
Matérn-3/2	1	$2.058 \times10^{-16}$	Fail	999	$5.091 \times10^{-16}$
Matérn-3/2	2	$2.040 \times10^{-16}$	Fail	998	$5.118 \times10^{-16}$
Matérn-3/2	5	$2.064 \times10^{-16}$	Fail	995	$5.099 \times10^{-16}$
Matérn-3/2	10	$2.073 \times10^{-16}$	Fail	990	$5.027 \times10^{-16}$
Matérn-5/2	0	$2.387 \times10^{-16}$	Pass	1000	$6.361 \times10^{-16}$
Matérn-5/2	1	$2.384 \times10^{-16}$	Pass	999	$6.239 \times10^{-16}$
Matérn-5/2	2	$2.386 \times10^{-16}$	Fail	998	$6.241 \times10^{-16}$
Matérn-5/2	5	$2.390 \times10^{-16}$	Fail	995	$6.258 \times10^{-16}$
Matérn-5/2	10	$2.403 \times10^{-16}$	Fail	990	$6.181 \times10^{-16}$
Matérn-7/2	0	$2.599 \times10^{-16}$	Pass	1000	$5.502 \times10^{-16}$
Matérn-7/2	1	$2.593 \times10^{-16}$	Fail	999	$5.489 \times10^{-16}$
Matérn-7/2	2	$2.643 \times10^{-16}$	Pass	998	$5.513 \times10^{-16}$
Matérn-7/2	5	$2.615 \times10^{-16}$	Fail	995	$5.466 \times10^{-16}$
Matérn-7/2	10	$2.607 \times10^{-16}$	Fail	990	$5.496 \times10^{-16}$

They are referred to as the polyharmonic and Matern operators, respectively. In addition to mathematical physics, these operators play an important role in statistical estimation using Gaussian random fields. A recurring computational kernel in these applications is the estimation of the Cholesky factors of their corresponding difference matrices.

However, discretizing these free-space PDEs requires us to make some assumptions about the boundary conditions of these operators [Wolhmuth]. Certain choices such as periodic or Neumann boundaries in certain domains can render the finite-difference matrix singular and in the case of certain higher-order operators (larger $m$ and $\nu$) numerically indefinite. We will demonstrate our method in this regime.

A second-order three-point stencil in one dimension is used to discretize the Laplacian in $\mathcal{L}_1,\mathcal{L}_2$ subject to periodic and Neumann boundary conditions. The relative reconstruction accuracies of matrices corresponding to different $m$ (correspondingly $\nu+ d/2$) using the finite-time methods and the pivoted Cholesky methods are compared in Table 8 and Table 9.

Operator	Boundary	Power	FT error	DPOTF2	Pivoted rank	Pivoted error
Polyharmonic	Neumann	1	$3.653 \times10^{-16}$	Fail	499	$1.526 \times10^{-15}$
Polyharmonic	Neumann	2	$4.277 \times10^{-16}$	Pass	499	$2.941 \times10^{-16}$
Polyharmonic	Neumann	3	$4.973 \times10^{-16}$	Fail	500	$1.744 \times10^{-16}$
Polyharmonic	Periodic	1	$3.793 \times10^{-16}$	Fail	499	$1.559 \times10^{-15}$
Polyharmonic	Periodic	2	$4.393 \times10^{-16}$	Pass	499	$4.272 \times10^{-16}$
Polyharmonic	Periodic	3	$4.524 \times10^{-16}$	Fail	500	$1.751 \times10^{-16}$
Matérn	Neumann	1	$4.056 \times10^{-15}$	Fail	499	$1.252 \times10^{-15}$
Matérn	Neumann	2	$4.053 \times10^{-15}$	Fail	499	$2.785 \times10^{-15}$
Matérn	Neumann	3	$5.756 \times10^{-15}$	Fail	499	$6.009 \times10^{-15}$
Matérn	Periodic	1	$3.978 \times10^{-15}$	Fail	499	$1.260 \times10^{-15}$
Matérn	Periodic	2	$4.188 \times10^{-15}$	Fail	499	$2.539 \times10^{-15}$
Matérn	Periodic	3	$5.682 \times10^{-15}$	Fail	499	$5.726 \times10^{-15}$

Taken together, these three sets of experiments strongly suggest that our finite-time method is a viable alternative to pivoted Cholesky implementations. A comprehensive HPC benchmark in the spirit of [Higham] to compare the latencies incurred due to pivoting and the associated gains made by the finite-time methods on modern computing architectures is an interesting direction for future work.

Discussion

This paper introduced the idea of finite-time regularization for factorizing positive semi-definite matrices. Canonical Cholesky factorization maps are not well defined in this regime. Therefore, the standard practice is to permute such matrices before applying this map. However, permutation incurs additional overhead, particularly for large matrices. Therefore, it is desirable to develop methods that don’t resort to pivoting. To this end, we propose finite-time Cholesky factorization.

The central idea is to design scalar dynamical systems whose finite-time evolution until $T$ provides a modified Cholesky map that is well-defined in the semi-definite case. The ambiguity of how to choose $T$ is clarified by deriving convergence rates for $\lVert A - L_T L_T^{\mathsf{T}} \rVert_F$ with respect to the evolution time $T$. Additionally, we have demonstrated with several numerical examples that finite-time methods are competitive with pivoted Cholesky methods. This provides an initial demonstration that the dynamical system interpretation provides a more stable Cholesky factorization method. Naturally, several question remain open.

Firstly, to facilitate a concise presentation, we focused on the properties of the proposed Cholesky map in infinite precision. Furthermore, all numerical experiments were carried out in dense double-precision floating-point representations of matrices. This leaves behind the obvious problem of establishing backwards stability of the finite-time Cholesky algorithm for a given floating-point representation. Naturally, one can also study how this algorithm behaves in mixed precision [13]. Understanding these questions will help us design more robust Cholesky methods in the future.

From a practical standpoint, when the matrices are sparse (as they often are), it is not well understood whether the finite-time Cholesky map will lower the number of pivots necessary to adequately factor that matrix. Furthermore, in cases where matrices have a low-rank structure, the finite-time Cholesky map can improve the robustness of Nyström methods [11] like RP-Cholesky [5], particularly when matrix representations could numerically lose positive-semi-definiteness. Empirical and theoretical investigations surrounding this can be consequential for several related problems in machine learning [14, 7, 6].

Lastly, establishing how the finite-time error and the numerical round-off errors induced during the computation of the Cholesky map; affect downstream applications [2] such as those where we solve linear equations, estimate the determinant of matrices, and similar matrix functions remains an active direction for future work.

Acknowledgments. RM was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.

Use of AI. The correctness of grammar in the text was verified with Grammarly. The numerical experiments were carried out initially by RM and ChatGPT Sol + Codex was used to wrap the initial research code into a working Julia library. The authors verified the correctness of the code in this library. The proofs for Theorems 2.1, 2.2 were originally derived by RM; ChatGPT Sol was used to simplify expressions, derive limits and polish the structure of the proof.

References