Interpolating log-determinant and trace of the powers of matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A} + t\textbf{B}$$\end{document}A+tB

We develop heuristic interpolation methods for the functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\mapsto \log \det \left( \textbf{A} + t\textbf{B} \right) $$\end{document}t↦logdetA+tB and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\mapsto {{\,\textrm{trace}\,}}\left( (\textbf{A} + t\textbf{B})^{p} \right) $$\end{document}t↦trace(A+tB)p where the matrices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}$$\end{document}A and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{B}$$\end{document}B are Hermitian and positive (semi) definite and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p$$\end{document}p and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t$$\end{document}t are real variables. These functions are featured in many applications in statistics, machine learning, and computational physics. The presented interpolation functions are based on the modification of sharp bounds for these functions. We demonstrate the accuracy and performance of the proposed method with numerical examples, namely, the marginal maximum likelihood estimation for Gaussian process regression and the estimation of the regularization parameter of ridge regression with the generalized cross-validation method.


Introduction
Estimation of the determinant and trace of matrices is a key component and often a computational challenge in many algorithms in data analysis, statistics, machine learning, computational physics, and computational biology. Some applications of trace estimation can be found in Ubaru and Saad (2018). A few examples of such applications are highperformance uncertainty quantification (Bekas et al. 2012;Kalantzis et al. 2013), optimal Bayesian experimental design (Chaloner and Verdinelli 1995), regression using Gaussian process (MacKay et al. 2003), rank estimation (Ubaru and Saad 2016), and computing observables in lattice quantum chromodynamics (Wu et al. 2016

Motivation
In this paper, we are interested in estimating the functions and where A and B are Hermitian and positive semi-definite (positive-definite if p < 0), and p and t are real numbers. 1 These functions are featured in a vast number of applications in statistics and machine learning. Often, in these applications, the goal is to optimize a problem for the parameter t, and the above functions should be evaluated for a wide range of t during the optimization process. A common example of such an application can be found in regularization techniques applied to inverse problems and supervised learning. For instance, in ridge regression by generalized cross-validation (Wahba 1977;Craven and Wahba 1978;Golub and von Matt 1997), the optimal regularization parameter t is sought by minimizing a function that involves (1b) at p = −1 (see Sect. 4.3). Another common usage of (1a) and (1b), for instance, is the mixed covariance functions of the form A + tI that appear frequently in Gaussian processes with additive noise (Ameli and Shadden 2022c, d) (see also Sect. 4.2). In most of these applications, the logdeterminant of the covariance matrix is common, particularly in likelihood functions or related variants. Namely, if one aims to maximize the likelihood by its derivative with respect to the parameter, the expression, ∂ ∂t log det(A + tI) = trace (A + tI) −1 , frequently appears. More generally, the function (1b) for p ∈ Z <0 appears in the | p|th derivative of such likelihood functions. Other examples of (1a) and (1b) are in experimental design (Haber et al. 2008), probabilistic principal component analysis (Bishop 2006, Sect. 12.2), relevance vector machines (Tipping 2001) and (Bishop 2006, Sect. 7.2), kernel smoothing (Rasmussen and Williams 2006, Sect. 2.6), and Bayesian linear models (Bishop 2006, Sect. 3.3).

Overview of related works
The difficulty of estimating (1a) and (1b) in all the above applications is that the matrices are generally large. Also, often in these applications, cases of particular interest in (1b) are when p < 0, but the | p|th inverse of the matrix A + tB is not available explicitly, rather it is implicitly known by matrix-vector multiplications through solving a linear system. Because of these, the evaluation of (1a) and (1b) are usually the main computational challenge in these problems, and several algorithms have been developed to address this problem.
The determinant and trace of the inverse of a Hermitian and positive-definite matrix can be calculated by the Cholesky factorization [cf. Eq. (24a) and (24b) in Sect. 4.2]. Using the Cholesky factorization, Takahashi et al. (1973) developed a method to find desired entries of a matrix inverse, such as its diagonals. The latter method was extended by Niessner and Reichert (1983). Also, Golub and Plemmons (1980) found entries of the inverse of the covariance matrix provided that the corresponding entries of its Cholesky factorization are non-zero. The complexity of this method is O(nw) where w is the bandwidth of the Cholesky matrix (see also Björck 1996, Sect. 6.7.4). Recently, probing and hierarchical probing methods were presented by Tang and Saad (2012) and Stathopoulos et al. (2013), respectively, to compute the diagonal entries of a matrix inverse.
In contrast to the above exact methods, many approximation methods have been developed. The stochastic trace estimator by Hutchinson (1990), which evolved from Girard (1989), uses Monte-Carlo sampling of random vectors with a Gaussian or Rademacher distribution. A similar concept was presented by Gibbs and MacKay (1997). Another randomized trace estimator was given by Avron and Toledo (2011) for symmetric and positive-definite implicit matrices. Based on the stochastic trace estimation, Wu et al. (2016) interpolated the diagonals of a matrix inverse. Also, Saibaba et al. (2017) improved the randomized estimation by a low-rank approximation of the matrix. Another tier of methods combines the idea of a stochastic trace estimator and Lanczos quadrature (Golub and Strakoš 1994;Bai et al. 1996;Bai and Golub 1997;Golub and Meurant 2009), which is known as stochastic Lanczos quadrature (SLQ). The numerical details of the SLQ method using either Lanczos tridiagonalization or Golub-Kahn bidiagonalization can be found for instance in Ubaru et al. (2017, Algorithms 1 and 2).

Objective and our contribution
Our objective is to develop a method to efficiently estimate (1a) or (1b) for a wide range of t. Note, if B is the identity matrix and the matrix A is small enough to pre-compute is eigenvalues, λ i (A), then, the evaluation of (1a) and (1b) for any t is immediate by However, for large matrices, estimating all eigenvalues is impractical. Thus, we herein develop an interpolation scheme for the functions (1a) and (1b) based on the following developments: • We present a Schatten-type operator that unifies the representation of (1a) and (1b) by a single continuous function. This operator leads to definitions of an associated norm and anti-norm on matrices. Sharp inequalities for this norm and anti-norm on the sum of two Hermitian and positive (semi) definite matrices provide rough estimates for (1a) and (1b). • We propose two interpolation methods based on the sharp norm and anti-norm inequalities mentioned above. Namely, we introduce interpolation functions based on a linear combination of orthogonal basis functions, or interpolation by rational polynomials.
We demonstrate the computational advantage of our method through two examples: • Gaussian process regression We compute (1a) and (1b) in the context of marginal likelihood estimation of Gaussian process regression. We show that with very few interpolation points, an accuracy of 0.01% can be achieved.
• Ridge regression We estimate the regularization parameter of ridge regression with the generalized crossvalidation method. We demonstrate that with only a few interpolation points, the ridge parameters can be estimated and the overall computational cost is reduced by 2 orders of magnitude.
The outline of the paper is as follows. In Sect. 2, we present matrix determinant and trace inequalities. In Sect. 3, we propose interpolation methods. In Sect. 4 we provide examples and a software package that implements the presented algorithms. In Sect. 5, we provide further applications of the method. Section 6 concludes the paper. Proofs are given in "Appendix A".

Determinant and trace inequalities
We will derive interpolations for (1a) and (1b) by modifying sharp bounds for these functions. In this section, we present these bounds. Without loss of generality, we temporarily omit the parameter t. However, in Sect. 3, we will retrieve the desired relations by replacing B with |t|B.
Let M n,m (C) denote the space of all n × m matrices with entries over the field C. We assume A, B ∈ M n,n (C) are Hermitian and positive semi-definite. Furthermore, for p < 0, we require matrices A and B to be positive-definite. The notations A B and A B on matrices A and B denotes A− B is positive-definite and positive semi-definite, respectively. Also, λ(A) := (λ 1 (A), . . . , λ n (A)) indicates the n-tuple of eigenvalues of matrix A.
Define a Schatten-class operator · p : M n,n (C) → R ≥0 by where |A| := √ A * A and A * denotes the Hermitian transpose of A. Since we assume the matrices are Hermitian and at least positive semi-definite, we omit |·| in subsequent expressions. Also, we note that p → · p is continuous at p = 0 since Namely, (4) is justified by observing that It is known that the generalized mean converges to the geometric mean, M 0 , as p → 0 (Hardy et al. 1952, p. 15), which concludes (4). For p ∈ [1, ∞), the operator · p is an equivalent norm to the Schatten p-norm of A. Conventionally, the Schatten norm is defined without the normalizing factor 1 n in (3), but the inclusion of this factor is justified by the continuity granted by (4). The Schatten norm is a subadditive function, meaning that it satisfies the triangle inequality The reverse triangle inequality follows from the above by provided that A B for (6b) to hold. For p = 1, the two above relations become equality by the additivity of the trace operator.
When p < 1, the operator (3) is no longer a norm, rather, is an anti-norm (Bourin and Hiai 2011) that satisfies superadditivity property A reversed inequality can also be derived from the above as provided that A B (or A B if p < 0) for (7b) to hold. We also note that inequality (7a) at p = 0 reduces to Brunn-Minkowski determinant inequality (Horn and Johnson 1990, p. 482, Theorem 7.8.8) For proofs and discussions of a general class of antinorms, which includes (3), we refer the reader to Hiai (2011, 2014). However, in "Appendix A", we provide a direct proof of (7a) and (7b) and the necessary and sufficient conditions for equality to hold in these relations for the operator (3) at p < 1.

Remark 1 (Comparisons to other inequalities)
There are other known bounds to functions (1a) and (1b). For instance, for the common case of p = −1, we can obtain the upper bound (Zhang 2011, p. 210, Theorem 7 Also, a lower bound can be obtained, for instance, by the arithmetic-harmonic mean inequality M −1 (λ(A ± B)) ≤ M 1 (λ(A ± B)), where M −1 and M 1 are the harmonic mean and arithmetic mean, respectively, (Mitrinović and Vasić 1970, Ch. 2, Theorem 1), which leads to The inequalities (9) or (10), however, are not as useful as the inequality (7a) for p = −1, since if B is either too small or too large compared to A, (9) and (10) do not asymptote to equality, whereas (7a) and (7b) become asymptotic equalities, which is a desired property for our purpose.

Interpolation of determinant and trace
We use the bounds provided by inequalities (6a), (6b), (7a), and (7b) to interpolate the functions (1b) and (1a). To this end, we replace the matrix B with |t|B in the bounds found in Sect. 2. Define We assume τ p,0 is known by directly computing ( 1 n trace (A p )) and (6b) and (7b) imply where t inf := inf{t |A + tB 0}. The above sharp inequalities become equality at t = 0. Also, (11a) and (11b) become asymptotic equalities as t → ∞. Based on the above, the bound function can be regarded as a reasonable approximation of τ p (t) at |t| τ p,0 where τ p (t) ≈ τ p,0 , and at t τ p,0 where τ p (t) ≈ t. We expectτ p (t) to deviate the most from τ p (t) when O(tτ −1 p,0 ) ≈ 1. Furthermore, to improve the approximation in the intermediate interval tτ −1 p,0 ∈ (c, c −1 ) for some c 1, we 2 Computing the determinant directly should be avoided as it can be a very large number. Rather, (det(·)) 1 n can be computed via exp( 1 n log det(·)). See Sect. 4.2 for an example.
define interpolating functions based on the above bounds to honor the known function values at some intermediate points t i ∈ (cτ p,0 , c −1 τ p,0 ). In particular, we specify interpolation points over logarithmically spaced intervals, because t is usually varied in a wide logarithmic range in most applications. We compute the function values at the interpolation points, τ p (t i ), with any of the trace estimation methods mentioned earlier.
Many types of interpolating functions can be employed to improve the above approximation. However, we seek interpolating functions whose parameters can be easily obtained by solving a linear system of equations. We define two such types of interpolations, namely, by a linear combination of basis functions and by rational polynomials, respectively in Sects. 3.1 and 3.2.

Interpolation with a linear combination of basis functions
Based on (12), we define an interpolating functionτ p (t) bỹ where φ i are basis functions and w i the weights. The basis functions for the domain t ∈ [0, ∞) can be used, which are inverse functions of the monomials and we refer to them as inversemonomials. These basis functions satisfy the conditions φ 0 (t) = t, φ i (0) = 0, and φ 0 (t) φ i (t), i > 0 when t 1. For consistency with (12), we set w 0 = 1. The coefficients w i , i = 1, . . . , q are found by solving a linear system of q equations using a priori known values τ p,i := τ p (t i ), i = 1, . . . , q. When q = 0, no intermediate interpolation point is introduced and the approximation function is the same as the boundτ p (t) given by (12).
Remark 2 An alternative could be to use monomials t i for interpolation functions, e.g., with w q+1 = 1, and the rest of the weights w i , i = 1, . . . , q determined from the known values of the function. This is not particularly useful in practice, as the exponentiation terms, t i , cause arithmetic underflows; also, Runge's phenomenon occurs even for low-order interpolations q > 1.
In practice, just a few interpolating points t i are sufficient to obtain a reasonable interpolation of τ p (t). However, when more interpolation points are used (such as when p ≥ 6), the linear system of equations for the weights w i becomes ill-conditioned. To overcome this issue, orthogonal basis functions can be used (see e.g., Seber and Lee 2012, Sect. 7.1 for a general discussion).
For our application, we seek basis functions φ ⊥ i (t) that are orthogonal on the unit interval t ∈ [0, 1]. Since we are interested in functions in the logarithmic scale of t, we define the inner product in the space of functions using the Haar measure d log(t) = dt/t . Applicability of Haar measure can be justified by letting t i = e x i , where x i are normally spaced interpolant points. Following the discussion of Seber and Lee (2012, Sect. 7.1) for linear regression using orthogonal polynomials, we use the conventional integrals with the Lebesgue measure dx to define the inner product of functions. The measure dx is equivalent to the Haar measure d log t for the variable t.
The desired orthogonality condition in the Hilbert space of functions on [0, 1] with respect to the Haar measure becomes where δ i j is the Kronecker delta function. A set of orthog- The first nine orthogonal basis functions are shown in Fig. 1 and the respective coefficients α i and a i j are given by Table  1. 3 A set of orthogonal functions can also be defined on intervals other than [0, 1] by adjusting the bounds of integration in (16), which yields a different set of function coefficients. However, it is more convenient to fix the domain of orthogonal functions to the unit interval [0, 1], and later scale the domain as desired, e.g., to [0, l] where l := max(t i ). Although this approach does not lead to orthogonal functions in [0, l], it nonetheless produces a well-conditioned system of equations for the weights w i .

Remark 3
The interpolation function defined in (13) asymptotes consistently toτ p (t) → t at t τ p,0 . On the other end, the convergenceτ p (t) → τ p,0 at t τ p,0 is not uniform, 3 We developed the python package ortho to generate an arbitrary number of orthogonal functions φ ⊥ j (t) using symbolic computations. See https://ameli.github.io/ortho for details.
Orthogonalized inverse-monomial functions in the logarithmic scale of t rather the interpolation function oscillates. This behavior is originated from the basis functions φ i , i > 0, that are not independent at t τ p,0 , particularly, near the origin. This dependency of basis functions cannot be resolved by the orthogonalized functions φ ⊥ i , as they are orthogonal with respect to the singular weight function t −1 at the origin. Thus, (13) should not be employed on very small logarithmic scales, rather, other interpolation functions should be employed for such purpose, such as presented in Sect. 3.2.

Interpolation with rational polynomials
We define another type of interpolating function that can perform well at small scales of t, by using rational polynomials. Definẽ which is the Padé approximation of τ p of order [q + 1, q]. We set a 0 = b 0 τ p,0 in order to satisfy τ p (0) = τ p,0 . Also, we note that the above interpolation satisfies the asymptotic relation τ p (t) → t as t → ∞. At q = 0, when no interpolant point is used, the above interpolation function falls back to (12) by setting b 0 = 1. For q > 0, 2q interpolant points t i are needed to solve the linear system of equations for the coefficients a 1 , . . . , a q and b 0 , . . . , b q−1 . An alternative rational polynomial is the Chebyshev rational function (Guo et al. 2002) where T i , i ∈ N are the Chebyshev polynomials of the first kind. The Chebyshev rational functions are orthogonal in [0, ∞) with respect to the weight function (t + 1) −1 √ t Chebyshev rational functions ). An interpolation of τ p (t) using Chebyshev rational functions can be given bỹ where α > 0 is a given scale parameter and will be explained shortly. Both sides of the above relation converge to zero at t → ∞. To satisfy τ The latter condition together with q linear equations on the interpolating points t i > 0, i = 1, . . . , q solve the weights w i . An advantage of using the above interpolation scheme is that we can arrange the interpolant points t i on the corresponding Chebyshev nodes to reduce the interpolation error. Figure 2 shows the Chebyshev rational basis functions in the form that are used on the right-hand side of (20). These basis functions converge at t 1 and t 1, whereas the main variability of these functions is mostly observed near O(t) ≈ 1. Thus, it is desirable to shift the interval of interpolation to the vicinity of O(t) = 1, which can be achieved by setting the scale parameter α. One approach to find an optimal interpolation parameter, α, is to minimize the curvature of the interpolating function, which is a common practice, for instance, in smoothing splines (Newbery and Garrett 1991). To this end, let w i (α) denote the solved weights for a given α. For simplicity, we transform the graph (t, Then, an optimal α * can be sought to minimize the arc integral of curvature squared of y α (x) by We note that in the absence of enough interpolant points, minimizing the curvature of the interpolating curve does not necessarily reduce interpolation error. However, when an adequate number of interpolant points are employed, the above approach can practically lead to a scale parameter α that enhances the interpolation. Finally, we note that unlike the interpolation scheme of Sect. 3.1 with the inverse monomial basis (14), both the Padé approximation of (18) and Chebyshev rational interpolation in (20) can interpolate τ p at negative values of t, namely in the domain t > t inf when t inf < 0 [see (11c) and (11d)].

Numerical examples
In Sect. 4.1, we briefly introduce a software package we developed for the presented numerical algorithm. This package was used to produce the results in Sects. 4.2 and 3.2. Indeed, the source code to reproduce the results and plots in the following sections can be found on the documentation of the software package. 4 Section 4.2 considers the problem of marginal likelihood estimation, which considers a full rank correlation matrix, and for this we use the interpolation functions of Sect. 3.1. Section 4.3 considers the problem of ridge regression, which considers a singular matrix, and for this we use the rational polynomial interpolation method of Sect. 3.2. We note that the interpolation with Chebyshev rational functions provide similar results to the orthogonalized inverse-monomials in (13) and we omit in our numerical examples for brevity.

Software package
The methods developed in this manuscript have been implemented into the python package imate, an implicit matrix trace estimator (Ameli and Shadden 2022b). This library estimates the determinant and trace of various functions of implicit matrices using either direct or stochastic estimation techniques and can process both dense matrices and large-scale sparse matrices. The main library of this package is written in C ++ and NVIDIA® CUDA and accelerated on both parallel CPU processors and CUDA-capable multi-GPU devices. The imate library is employed in the python package glearn, a machine learning library using Gaussian process regression (Ameli and Shadden 2022a).
Interpolant points Interpolated Fig. 3 Result of the code in Lisiting 1. a Comparison of interpolated versus exact value of the function τ −1 (t). The exact function (red curve) is overlaid by the interpolated curve (black curve). b Relative error of the comparison lates f p : t → A + tB p . Briefly, Line 9 generates a sample correlation matrix A ∈ M n,n (R) on a randomly generated set of n = 50 2 points using an exponential decay kernel. In Line 15, we create an instance of the class imate.InterpolateSchatten. Setting B=None indicates B is the identity matrix using an efficient implementation that does not require storing identity matrix. The instantiation in Line 19 internally computes τ −1,i = τ −1 (t i ) on eight interpolant points t i = 10 −4 , 10 −3 , . . . , 10 3 and obtains the interpolation coefficients for the orthogonalized inverse-monomial basis functions (14) since kind='IMBF' was specified. Other possible methods can be the exact method with no interpolation (EXT), eigenvalue method (EIG) given in (2b), monomial basis functions (MBF) given in (15), Padé rational polynomial functions (RPF) given in (18), Chebyshev rational functions (CRF) given in (20), or radial basis functions (RBF) (which we do not cover herein for brevity). The evaluation of τ −1,i can be configured by passing a dictionary of settings to the options argument, and we refer the interested reader to the package documentation for further details. In this minimalistic example, we compute τ −1,i using Cholesky decomposition, as further detailed in Sect. 4.2. Other methods include stochastic Lanczos quadrature or Hutchinson estimation; we compare such methods in Sect. 4.3. Once the interpolation object is initialized, future calls to interpolate an arbitrary number of points t are returned almost instantly. In Line 19, the interpolation is performed on 1000 points in the interval [10 −4 , 10 3 ] spaced uniformly on the logarithmic scale. A comparison of the interpolated result versus the exact solution are shown in Fig. 3. It can be seen that with only eight interpolant points, the relative error of interpolation over a wide range of parameter t is around 0.1%.

Marginal likelihood estimation for Gaussian process regression
Here we generate a full rank correlation matrix from a spatially correlated set of points x ∈ D = [0, 1] 2 . To define a spatial correlation function, we use the isotropic exponential decay kernel given by where ρ is the correlation scale, set to ρ = 0.1. The above exponential decay kernel represents an Ornstein-Uhlenbeck random process, which is a Gaussian and zeroth-order Markov process (Rasmussen and Williams 2006, p. 85). To produce discrete data, we sample n = 50 2 points from D, which yields the symmetric and positive-definite correlation matrix A with the components A i j = κ(x i , x j |ρ). We aim to interpolate functions for p = −1, −2, which appear in many statistical applications, such as the estimation of noise in Gaussian process regression (Ameli and Shadden 2022c). Specifically, the above functions for p = 0, −1, and −2 appear in the corresponding likelihood function, and its Jacobian and Hessian, respectively. We compute the exact value of τ p (t) for p ∈ Z ≤0 (either at interpolant points t i or at all points t for the purpose of benchmark comparison) as follows. We compute the Cholesky factorization of (A + tI) | p| = L | p| L | p| , where L | p| is lower triangular. Then where (L 1 ) ii is the ith diagonal element of L 1 and · F is the Frobenius norm. The second equality in (24b) employs the cyclic property of the trace operator. A simple method to compute L −1 | p| 2 F without storing L −1 | p| is to solve the lower triangular system L | p| x i = e i for x i , i = 1, . . . , n, where e i = (0, . . . , 0, 1, 0, . . . , 0) is a column vector of zeros, except, its ith entry is one. The solution vector x i is the ith This method is memory efficient since the vectors x i do not need to be stored.
We note that the complexity of the interpolation method is the number of evaluations of τ p at interpolant points t i and at t = 0 (which is proportional to q) times the complexity of computing τ p at a single point t. For instance, by using the Cholesky method in (24a) or (24b) which costs O( 1 3 n 3 ) for a matrix of size n, the complexity of the interpolation method is O( 1 3 n 3 q). Remark 4 (Case of Sparse Matrices) There exist efficient methods to compute the Cholesky factorization of sparse matrices (see e.g., Davis 2006, Ch. 4). Also, the inverse of the sparse triangular matrix L | p| can be computed at O(n 2 ) complexity (Stewart 1998, pp. 93-95), and a linear system with both sparse kernel L | p| and sparse right-hand side e i can be solved efficiently (see Davis 2006, Sect. 3.2).
The exact value of τ p (t), for p = 0, −1, −2, computed directly using the Cholesky factorization method described above are respectively shown in Fig. 4a, c, e by the solid black curve (overlaid by the red curve) in the range t ∈ [10 −4 , 10 3 ]. The dashed black curves in Fig. 4a, c, e are the lower boundŝ τ p (t) given by (12), which can be thought of as the estimation with zero interpolant points, i.e., q = 0. For completeness, we have also shown the upper bound of τ −1 (t) by the black dash-dot line in Fig. 4c, given by The above upper bound can be obtained from (10) and the fact that trace(A) = n, since the diagonals of the correlation matrix are 1. However, unlike the lower bound in (7a), the upper bound (25) is not useful for approximation as it does not asymptote to τ −1 (t) at small t. Nonetheless, both the lower and upper bounds asymptote to t at large t.
To estimate τ p , we used the interpolation function in (13) with the set of orthonormal basis functions in Table 1. The colored solid lines in Fig. 4a, c, e are the interpolationsτ p (t) with q = 1, 3, 5, 7, and 9 interpolant points, t i , spanning from 10 −4 to 10 3 . It can be seen from the embedded diagrams in Fig. 4a, c, e thatτ p (t) is remarkably close to the true function value. In practice, fewer interpolant points in a small range, e.g., [10 −2 , 10 2 ], are sufficient to effectively interpolate τ p .
To better compare the exact and interpolated functions, the relative error of the interpolations is shown in Fig. 4b, d, f. The relative error of the lower bound (dashed curve) rapidly vanishes at both ends, namely, at t τ p,0 and t τ p,0 , where τ 0,0 = 0.22, τ −1,0 = 0.16, and τ −2,0 = 0.14. The absolute error of the upper bound is highest at O(tτ −1 p,0 ) = 1, or t ≈ τ p,0 , which is slightly to the left of the relative error peak on each diagram.
Based on the lower bound, we distribute the interpolant points, t i , almost evenly around t ≈ τ p,0 where the lower bound has the highest error. The blue curve in Fig. 4b, d, f corresponds to the case with only one interpolation point at t 1 = 10 −1 , which already leads to a relative error less than 3% almost everywhere. On the other hand, with only nine interpolation points t i ∈ {10 −4 , 4 × 10 −4 , 10 −3 , 10 −2 , . . . , 10 3 } the relative error becomes less than 0.01%. Beyond the strong The interpolations using 7 and 9 points lead to relative errors of less than 0.02% and 0.01%, respectively. Rows correspond to p = 0, −1, and −2, respectively accuracy shown by the relative errors, the absolute errors are more compelling since τ p (t) decays by orders of magnitude at large t, making the absolute error negligible at t τ p,0 .

Ridge regression with generalized cross-validation
Here we calculate the optimal regularization parameter for a linear ridge regression model using generalized crossvalidation (GCV). Consider the linear model y = Xβ + , where y ∈ R n is a column vector of given data, X ∈ M n,m (R) is the known design matrix representing m basis functions where m < n, β ∈ R m is the unknown coefficients of the linear model, and ∼ N (0, σ 2 K) is the correlated residual error of the model, which is a zero-mean Gaussian random vector with the symmetric and positivedefinite correlation matrix K and unknown variance σ 2 . A generalized least-squares solution to this problem minimizes the square Mahalanobis distance y − Xβ 2 K −1 := ( y − Xβ) K −1 ( y − Xβ) yielding an estimation of β bŷ β = (X K −1 X) −1 X K −1 y (Seber and Lee 2012, p. 67).
When X is not full rank, the least-squares problem is not well-conditioned. A resolution of the ill-conditioned problems is the ridge (Tikhonov) regularization, where the function y − Xβ 2 K −1 + nθ β 2 is minimized instead (Seber and Lee 2012, Sect. 12.5.2). Here, the penalty term is β 2 = β β where is the symmetric and positivedefinite penalty matrix. The estimate of β using the penalty term becomeŝ Also, the fitted values on the training points areŷ = Xβ, which can be written asŷ = S θ y, where the smoother matrix S θ is defined by The regularization parameter, θ , plays a crucial role to balance the residual error versus the added penalty term. The generalized cross-validation method (Wahba 1977;Craven and Wahba 1978;Golub et al. 1979) is a popular way to seek an optimal regularization parameter without needing to estimate the error variance σ 2 . Namely, the regularization parameter is sought as the minimizer of (Hastie et al. 2001, p. 244). 5 For large matrices, it is difficult is to compute trace(S θ ) (also known as the effective degrees of freedom) in the denominator of (28), and several methods have been developed to address this problem (Golub and von Matt 1997;Lukas et al. 2010).

Estimating the trace
Using the presented interpolation method, we aim to estimate trace(S θ ). Let = LL be the Cholesky decomposition of . Using the cyclic property of trace operator, we have In the above, I m×m is identity matrix of size m. To compute the above term, we interpolate where t := nθ − s and We note that the size of A and I is m. Also, A is symmetric and positive-definite since it can be written as a Gramian matrix. The purpose of the fixed parameter s 1 is to slightly shift the singular matrix L −1 X K −1 XL − to make A non-singular. The shift is necessary since without it, (30) is undefined at t = 0, and we cannot compute τ −1,0 = m/ trace(A −1 ). Also, the shift can improve interpolation by relocating the origin of t to the vicinity of the interval where we are interested to compute V (θ ).
For simplicity in our numerical experiment, we set K and to identity matrices of sizes n and m, respectively. We also set s = 10 −3 . We create an ill-conditioned design matrix X for our numerical example using singular value decomposition X = U V . The orthogonal matrices U ∈ M n,n (R) and V ∈ M m,m (R) were produced by the Householder matrices where u ∈ R n and v ∈ R m are random vectors (see also Golub and von Matt 1997, Sect. 10). The diagonal matrix ∈ M n,m (R) was defined by We set n = 10 3 and m = 500. We generated data by letting y = Xβ + , where β, and , were randomly generated with a unit variance, and σ = 0.4, respectively.
We computed the exact solution of τ −1 (t) in (30), and interpolation points τ −1 (t i ), using the Cholesky factorization method described by (24b). The exact solution is shown by the solid black curve in Fig. 5 (overlaid by the red curve) with τ −1,0 = 960.5 −1 . The lower boundτ −1 (t) from (7a) is shown by the dashed black curve for t > 0. In contrast, at t ∈ (t inf , 0], the upper bound from (7b) is shown, where t inf = − min(λ(A)) and min(λ(A)) ≈ s = 10 −3 is the smallest eigenvalue of A. The relative error of the bounds with respect to the exact solution are shown in Fig. 5b. The peak of the absolute error of the lower bound is located approximately at t ≈ τ −1,0 ≈ 10 −3 , and the peak of its relative error of the lower bound is slightly to the right of this value.
The interpolation functionτ −1 (t) with q = 1, 2, 3 is shown in Fig. 5a using 2q interpolation points t i in the interval t i ∈ [5 × 10 −3 , 5] that are equally distanced in the logarithmic scale. The red curve corresponding to q = 3 and the black curve (exact solution) are indistinguishable even in the embedded diagram that magnifies the location with the highest error. The relative error of the interpolations is shown by Fig. 5b. On the far left of the range of t, the error spikes due to the singularity of the matrix X, which makes τ −1 (t) undefined at t = −10 −3 , corresponding to θ = 0. On the rest of the range, the green and red curves respectively show less than 0.1% and 0.05% relative errors, which are compelling accuracy for a broad range of t, and achieved with only four and six interpolation points, respectively.

Optimization of generalized cross-validation
Here we apply the result of our trace interpolations above to solve the generalized cross-validation problem. The function V (θ ) from (28) is plotted in Fig. 6, with the black curve, corresponding to the exact solution with τ −1 (t) applied in the denominator of V (θ ), serving as a benchmark for comparison. The blue, green, and red curves correspond to the proceeding trace interpolations applied in the denominator of V (θ ). The interpolated curves exhibit both local minima of V (θ ) similar to the benchmark curve, but with slight differences in the positions of the minima. Due to the singularity at θ = 0, the interpolations of τ −1 (t) become less accurate

Generalized Cross Validation
No interpolation (Exact) Interpolation, q = 1 Interpolation, q = 2 Interpolation, q = 3 Fig. 6 The generalized cross-validation function is shown, where the black and colored curves correspond respectively to the exact and interpolated computation of τ −1 (t) in the denominator of V (θ). The global minimum of each curve is shown by a dot at low values of θ . At higher values of θ , all curves steadily asymptote to a constant. We note that the results in Fig. 6 are compelling since the estimation of V (θ ) is sensitive to the interpolation of its denominator. Namely, a consistent interpolation accuracy over all the parameter range is essential to capture the qualitative shape of V (θ ) correctly.
The global minimum of V (θ ) at θ = θ * is the optimal compromise between an ill-conditioned regression problem (small θ) and a highly regularized regression problem (large θ ). We aimed to test the practicality of our interpolation method in searching the global minimum, V (θ * ), by numerical optimization. We note that we constructed X in (31) so that V (θ ) would have two local minima thus making optimization less trivial. In general, the generalized cross-validation function may have more than one local minimum necessi-tating global search algorithms (Kent and Mohammadzadeh 2000). The optimization was performed using a differential evolution optimization method (Storn and Price 1997) with a best/1/exp strategy and 40 initial guess points. The results are shown in the first four rows of Table 2, where the trace of a matrix inverse is computed by the Cholesky factorization described in (24b). In the first row, τ −1 is computed exactly, i.e., without interpolation, at all requested locations t during the optimization procedure. On the second to fourth rows, τ −1 is first pre-computed at the interpolation points, t i , by the Cholesky factorization, and then the interpolation is subsequently used during the optimization procedure.
In the table, N tr counts the number of exact evaluations of τ −1 . For the Padé rational polynomial interpolation method of degree q, we had N tr = 2q + 1, accounting for 2q interpolant points in addition to the evaluation of τ −1,0 at t = 0. N tot is the total number of estimations of τ −1 during the optimization process. In the first row, N tr = N tot as all points are evaluated exactly, i.e., without interpolation. However, for the interpolation methods, N tot consists of N tr plus the evaluations of τ −1 via interpolation.
The exact computations of τ −1 (at N tr points) are the most computationally expensive part of the overall process. Our numerical experiments were performed on the Intel Xeon E5-2640 v4 processor using shared memory parallelism. We measured computational costs by the total CPU processing time of all computing cores. T tr denotes the processing time of computing τ −1 exactly at the N tr points. T tot measures the processing time of the overall optimization, which includes T tr . As shown, the interpolation methods took significantly less processing time compared to no interpolation, namely, by two orders of magnitude for T tr , and an order of magnitude for T tot . We also observe that without interpolation, T tr is the Table 2 Comparison of methods to optimize the regularization parameter θ, with and without interpolation of τ −1 (t), and by various algorithms of computing trace of a matrix inverse dominant part of the total processing time, T tot . In contrast, with interpolation, T tr becomes so small that T tot is dominated by the cost of evaluating the numerator of V (θ ) in (28), which is proportional to N tot . The results of computing the optimized parameter, θ * , and the corresponding minimum, V (θ * ), are shown in the seventh and eighth columns of Table 2, respectively. The ninth column is the relative error of estimating θ * , and obtained by comparing log 10 θ * between the interpolated and benchmark solution (i.e., first row). The last two columns are the 2 norm of the error ofβ [using (26)] andŷ compared to their exact solution, and their relative error are obtained by normalizing with the 2 norm of their exact solution. We observed, for example for q = 2, that with one-tenth of the processing time, an accuracy of 2% error forŷ is achieved, which is generally sufficient in practical applications. Also, for q = 3, the error reduces to <1% with similar processing time. In general, the error can be improved simply by using more interpolant points. We have found that simple heuristics for setting defaults for q are broadly effective. Namely, if θ * is expected to be found in an known interval, one can use a small number of interpolating points (q = 1 ∼ 2) on the boundary or center of the interval. If there is no prior knowledge of the range, one can let an optimization scheme search for θ * in a wide logarithmic range, e.g., [10 −7 , 10 +1 ] with q = 3 ∼ 4.

Testing alternative trace estimators
Besides the Cholesky factorization algorithm, we also repeated the numerical experiments above with stochastic trace estimators, namely, the Hutchinson's algorithm (Hutchinson 1990) and stochastic Lanczos quadrature algorithm (Golub and Meurant 2009, Sect. 11.6.1) to compute the trace of a matrix inverse. These class of randomized methods are attractive due to their scalability to very large matrices, where employing the exact methods could be inefficient, if not infeasible. However, these methods do not compute the exact value of determinant or trace, rather they converge to a solution by Monte-Carlo sampling through iterations. The complexity of Hutchinson method using conjugate gradient is O(ρn 2 s) where ρ is the density of matrix (ρ = 1 for dense matrices) and s is the number of random vectors for Monte-Carlo sampling. We recall that in our application, the cost of the interpolation scheme is q times the above-mentioned complexity, i.e., O qρn 2 s .
Alternatively, the computational cost of the SLQ method is O((ρn 2 + nl)sl) where l is the Lanczos degree, which is the number of Lanczos tri-diagonalization iterations (see details e.g., in Ubaru et al. 2017, Sect. 3). Thus, the complexity of the interpolation method becomes O q(ρn 2 + nl)sl .
In both these algorithms, we employed s = 30 random vectors with Rademacher distribution for Monte-Carlo sampling. Also, in SLQ algorithm, we set the Lanczos degree to 30.
The results for Hutchinson's algorithm are shown in the fifth to eighth rows, and results for the SLQ algorithm are shown in the ninth to twelfth rows, of Table 2. Similar to the Cholesky factorization, in both stochastic estimators, the interpolation technique reduces the processing times compared to no interpolation, namely, T tr is reduced by two orders of magnitude, and T tot by an order of magnitude, while maintaining a reasonable accuracy.
We note that the interpolation with the stochastic methods introduces error due to the uncertainly in the randomized estimation of τ −1 at the interpolant points t i . However, the additional error caused by interpolation itself can be less than the error due to aforementioned stochastic estimation. For instance, without interpolation, the SLQ method estimateŝ y with a 1.58% error, whereas, interpolation with q = 3 results in a similar error of 2.76% but at a greater than 20fold reduction in computational cost.

Further applications
We recall that the presented interpolation scheme can be applied to any formulation that consists of the trace or determinant of a power of the one-parameter affine matrix function A + tB where A and B are Hermitian and positive-definite. Often in applications, an algebraic trick [such as in (29)] is required to form such an affine matrix function. We here provide two other closely related examples where such affine matrix function can be formulated.

Reproducing kernel Hilbert space
Let H K be a reproducing kernel Hilbert space equipped with the reproducing kernel K that defines the function evaluation f (x) = f , K (·, x) H K . Consider an infinite-dimensional generalized ridge regression on H K to estimate y = f (x) with the given training set {(x i , y i )} n i=1 by the minimization problem (Hastie et al. 2001, Sect. 5.8.2) The solution to the above problem has the form f (·) = j α j K (·, x j ). For the finite-dimensional formulation, define the kernel matrix K with the components K i j := K (x i , x j ), which is symmetric and positive-definite. Let α := [α 1 , . . . , α n ] and y := [y 1 , . . . , y n ] . The minimization problem in finitedimensional setting becomes where α 2 K = α Kα. The optimal solution to the above problem iŝ α = (K + θ I) −1 y, and the fitted values on the training points areŷ = Kα =: S θ y where the smoother matrix S θ is defined by S θ := K(K+ θ I) −1 .
One may seek the optimal value for θ as the minimizer of the GCV function We recall that the expensive part of computing (32) is the term trace(S θ ). To apply our interpolation scheme, write S θ as the Reinsch form We realize that where t := θ −1 . The proposed interpolation method follows by using A = K −1 , B = I, and τ −1,0 = n/ trace(K).

Kernel-based GCV for mixed models
Another formulation of kernel-based GCV, for instance by Xu and Zhu (2009, Eqs. 9 and 10), yields a function of the form In the above, the covariance = (φ) is symmetric and positive-definite, the design matrix of random effects Z has full column-rank, andH =H(h) is the smoother matrix when the random effects are absent. Optimal values of the parameters (h, φ) are sought by minimizing V .
It is possible to represent the term in the denominator of (33) by the trace of the inverse of a single matrix to be written as τ −1 . To do so, let P := I −H and Y := PZ. Using the Woodbury matrix identity and (34), we can represent the term inside the trace in (33) as If P is positive-definite, then let P = LL be the Cholesky decomposition of P. By using the cyclic property of trace operator, we have Note that both matrices A := L −1 L − and B := L −1 Y −1 Y L − are symmetric and positive-definite since they are in the Gramian form. To compute (35), the presented interpolation method can be applied for instance if (φ) is linear in its parameter. Such assumption is common, for instance when = φK where φ is variance and K is the correlation matrix. In such a case, the sum of two matrices in (35) becomes an affine function of t := φ −1 and trace(I −H(h, φ)) can be written as τ −1 (t).

Conclusions
In many applications in statistics and machine learning, it is desirable to estimate the determinant and trace of the real powers of a one-parameter family of matrix functions A+tB where the parameter t varies and the matrices A and B in the formulation remain unchanged. There exist many efficient techniques to estimate the determinant and trace of implicit matrices (such as inverse of a matrix), however, these methods are geared toward generic matrices. Using those methods, the computation of the determinant and trace of the parametric matrices should be repeated for each parameter value as the matrix is updated. To efficiently perform such computation for a wide range of parameter t, we presented in this work heuristic methods to interpolate the functions t → log det(A + tB) and t → trace((A + tB) p ). The interpolation approach is based on sharp bounds for these functions using inequalities for a Schatten-type norm and anti-norm. We proposed two types of interpolation functions, namely, interpolation with a linear combination of orthogonalized inverse-monimial basis functions, and interpolation with rational polynomials, which includes Padé approximation and Chebushev rational functions. We demonstrated that both functions can provide highly accurate interpolation over a wide range of t using very few interpolation points. The rational polynomials generally provide better results in the neighborhood of the origin of the parameter. In the regions away from the origin, choice of interpolation method is less important; namely we observed e.g., the interpolation with Chebyshev rational functions provide similar results to the orthogonalized inverse-monomials in (13) in such cases. All the presented interpolation methods can lead to one to two orders of magnitude savings in processing time in practical applications that require frequent evaluations of log det(A + tB) or trace((A + tB) p ).
For applications where one is interested in values of t τ p,0 (such as in Sect. 4.3 where the matrix was shifted due to being ill-conditioned) interpolation using (18) is recommended. One should keep in mind that there exists the possibility that (18) can become singular at its poles, but a slight rearrangement of the interpolant points t i can be used to ensure these poles are outside the domain of interest. Although (18) provides accurate interpolation for a broad range of t, for a higher number of interpolation points (e.g 6 or more), relation (13) or (20) is preferred.
In closing, the presented interpolation method can be effectively utilized on large data, particularly with the powerful framework of randomized estimators of trace and log-determinant. A practical application of this method together with stochastic Lanczos quadrature on sparse matrices is given by Ameli and Shadden (2022c) to efficiently train Gaussian process regression. The interested reader may refer to Ameli and Shadden (2022b) where the interpolation scheme can be applied to massive data (e.g., n ∼ 2 25 ) using the imate package. We prove Theorem 1 as follows.
Definition 1 (Majorization) For the n-tuple x = (x 1 , . . . , x n ) ∈ R n , we denote by x ↓ = (x ↓ 1 , . . . , x ↓ n ) the tuple with the same components as x, but sorted in decreasing order. Let x, y ∈ R n . We say x weakly majorizes y from below (submajorizes) and indicate by x w y if and only if Proof We proceed the proof for p < 0 as the case p ∈ (0, 1) follows similarly. By Ky Fan eigenvalue inequality for Hermitian matrices (Zhang 2011, p. 356 the equality in Proposition 2 holds if A and B commute, which means A and B have the same eigenspace (Horn and Johnson 1990, p. 50, Theorem 1.3.12). By combining these two conditions, equality in (A.1a) and (A.1b) is achieved when A is a scalar multiple of B.