OPTIMIZATION OVER COVARIANCE MATRICES WITH A PARAMETERIZED METRIC Yibang Li1 , Bamdev Mishra2 , Pratik Jawanpuria3 and Cyrus Mostajeran1 1
arXiv:2609.17089v1 [math.OC] 15 Sep 2026
2
Nanyang Technological University, Singapore 3 Microsoft India Indian Institute of Technology Bombay, India [email protected] [email protected] [email protected] [email protected] ABSTRACT
The choice of Riemannian metric can strongly influence the convergence of gradient-based optimization over covariance matrices. Euclidean, Bures–Wasserstein and affine-invariant metrics are common choices, but their relative effectiveness depends on the objective. We introduce a two-parameter family defined by X p LX q + X q LX p = U , solved for L at each tangent vector U , that contains all three as exact members, at (0, 0), (1, 0) and (1, 1), and extends past them. We treat the choice of member as a particular way of preconditioning for a given problem. To this end, we analyze the conditioning of the Riemannian Hessian at the solution. We show that it obeys a lower bound that depends on (p, q) only through the exponent r = p + q. When the Euclidean Hessian is a pure power that mixes no eigendirections, the member p = q = r/2 attains that bound, and a closed-form criterion identifies the other members that do. We discuss ways to tune r for a given problem. Experiments on real covariance data confirm the predicted conditioning and the benefit of tuning r. A task covariance example shows a further gain from tuning the shape. Index Terms— Covariance matrices, Riemannian optimization, preconditioning, Bures–Wasserstein, metric learning
[12] who compare Bures–Wasserstein (BW) and affine-invariant geometry and report that each wins on a different objective class, while [13] introduce a parameterized generalization. More concretely, preconditioning on manifolds has attracted much attention. Mishra and Sepulchre [11] build the metric from the Hessian of the Lagrangian. Shustin and Avron [14] choose the metric on the generalized Stiefel manifold through a preconditioning scheme, motivate the choice by the condition number of the Riemannian Hessian at the optimum, and identify the ideal preconditioner as the Euclidean Hessian there. Gao et al. [15] do the same on product manifolds. Closest to this work, Zhou et al. [16] discuss Alpha-Procrustes metrics (which include logEuclidean and BW) for optimization, especially from a robustness viewpoint. Contributions. We explore the question of metric choice for the SPD manifold. To this end, our contributions are the following. We introduce a two-parameter metric family containing Euclidean, Bures–Wasserstein and affine-invariant geometry as exact members. We show how the condition number of the Hessian depends on p and q, and motivate ways to approximate this. Finally, our experiments show the benefit of tuning the metric. 2. A PARAMETERIZED METRIC
1. INTRODUCTION Symmetric positive-definite (SPD) matrices are the working object in covariance estimation, adaptive beamforming, radar detection, diffusion tensor imaging [1], and brain–computer interfaces [2]. The set of SPD matrices forms a manifold Sn ++ , an open subset of the space Sn of symmetric matrices, so Sn is its tangent space at every point. Endowed with a metric, it has the structure of a Riemannian manifold. Minimizing a function f over Sn ++ by gradient descent requires a choice of Riemannian metric, which determines the gradient and hence the descent direction. Each metric in common use on the set of SPD matrices is often argued for on geometric grounds. Affine-invariant geometry is complete and congruence invariant. Bures–Wasserstein geometry is the covariance part of quadratic optimal transport between Gaussian measures [3, 4, 5]. Log-Euclidean geometry is widely used in imaging [6]. Thanwerdas and Pennec [7] classify the O(n)-invariant metrics and show that the kernel metrics of Hiai and Petz [8], each fixed by a single function of two eigenvalues, form a subclass containing most of the metrics in common use, and elsewhere they build continua that interpolate between named metrics [9, 10]. Since the metric is what converts ∇f into a search direction, it acts as an intrinsic preconditioner rather than a passive modeling choice [11]. Picking one is therefore an algorithmic decision rather than a geometric one. This is, for example, explored by Han et al.
The metric. Fix two real numbers p and q, and let X ∈ Sn ++ be the point at which the metric is being defined. This follows the usual definition of Bures–Wasserstein geometry [5]. For a tangent vector U ∈ Sn let L = Lp,q (U ) solve the Sylvester-type equation X p LX q + X q LX p = U.
(1)
For tangent vectors U and V , the Riemannian metric g is defined as 2 (p,q) gX (U, V ) = 21 cp,q tr Lp,q (U ) V , cp,q = 4 1−(p−q) , (2) where cp,q is a normalizing constant. Let d1 , . . . , dn and P hold the eigenvalues and eigenvectors of X, and write M ′ = P ⊤ M P for the eigenbasis coordinates of any M ∈ Sn . In this eigenbasis, the operator on the left-hand side of (1) acts entrywise, multiplying L′ij by dpi dqj + dqi dpj . Since the eigenvalues of X are positive, these coefficients are strictly positive for all real p and q. The operator is therefore self-adjoint, positive definite, and invertible, so (2) defines a Riemannian metric for every such pair. At (p, q) = (1, 0), (1) reduces to XL + LX = U , and (2) recovers the Bures–Wasserstein metric. The weight. Solving (1) entry by entry puts (2) in the kernel form of Hiai and Petz [8], ′ X Uij Vij′ ϕ gX (U, V ) = , ϕ = ϕp,q , (3) ϕ(di , dj ) i,j
Algorithm 1 Descent on Sn ++ in the metric (2)
with the weight in the denominator, 2 ϕp,q (x, y) = 12 4 (p−q) xp y q + xq y p ,
(4)
symmetric in (p, q) and positively homogeneous of degree r = p+q, meaning ϕp,q (tx, ty) = t r ϕp,q (x, y) for every t > 0. We call r the exponent of the member of the metric family. Named members. Evaluating (4) at (p, q) = (0, 0), (1, 0) and (1, 1) gives the weights 1, 2(x + y) and xy, of degrees 0, 1 and 2. These are exactly the Euclidean, Bures–Wasserstein and affineinvariant metrics, respectively. The metric family interpolates between them, and since p and q may be any reals, it extends past them in every direction. On the diagonal p = q = r/2 (1) reads 2X r/2 LX r/2 = U , so L = 12 X −r/2 U X −r/2 and (2) becomes the power family in closed form, ϕ gX (U, V ) = tr X −r/2 U X −r/2 V , ϕ(x, y) = (xy)r/2 . (5) Relation to known metric families. Every smooth positive ϕ defines a kernel metric [8], so (4) is a subfamily of a known class, singled out by the factorization ϕp,q (x, y) = 4
(p−q)2
(xy)
r/2
cosh p−q log(x/y) 2
,
(6)
in which r fixes the homogeneity degree and p − q the shape at fixed degree. That separation is what makes the conditioning law of Section 3 possible. It is the separation rather than any individual member that is new. The family also shares members with the mixedpower-Euclidean family of [10]. The inverse-Euclidean metric is the pullback of the Euclidean metric under X 7→ X −1 . It is our diagonal member at r = 4 and the mixed-power-Euclidean member at (−1, −1). Off the diagonal, our member (p, q) = ( 12 , 0) agrees up to a constant factor with their member at (1, 12 ). The family differs from the generalized Bures–Wasserstein geometry of [13], which deforms the metric by an SPD matrix parameter rather than by two scalars. Descent on the metric family. Write G for the metric operator ϕ defined by gX (U, V ) = tr(G[U ]V ). By (3) it divides entrywise in the eigenbasis, so its inverse multiplies entrywise, G −1 [Z] = P (K
Z ′ )P ⊤ ,
Kij = ϕp,q (di , dj ),
(7)
where is the entrywise product. The Riemannian gradient is grad f = G −1 [G] = P (K G′ )P ⊤ for G = ∇f (X). The member (p, q) enters Algorithm 1 only through the line that forms K, so one implementation covers the whole family at the cost of the eigendecomposition already incurred. For the objectives below that order is already paid to form G, so the metric adds no extra cost. Where the gradient is cheaper, the eigendecomposition sets the cost and limits the reachable n. Steps are taken with the retraction RX (V ) = X + V + 12 V X −1 V,
(8)
which stays positive definite for every symmetric V . It is a retraction in the usual sense, since RX (0) = X and DRX (0)[V ] = V . Most members have no exponential map in closed form, and (8) replaces it at the price of one triple product. 3. THE ROLE OF r = p + q IN HESSIAN CONDITIONING Finding the right (p, q) for a given problem looks like a twodimensional search, scored by the conditioning of the Riemannian Hessian of the objective in the new metric. Below, we give a principled way to choose r, p and q.
Require: f , X0 ∈ Sn ++ , member (p, q), steps tk 1: for k = 0, 1, 2, . . . do 2: G ← ∇f (Xk ) ▷ Euclidean gradient 3: (d, P ) ← eig(Xk ) ▷ O(n3 ) 4: Kij ← ϕp,q (di , dj ) ▷ weight (4) 5: U ← P K (P ⊤ GP ) P ⊤ ▷ Riemannian gradient (7) t2
6: Xk+1 ← Xk − tk U + 2k U Xk−1 U 7: end for
▷ retract (8)
The useful quantity is the condition number κϕ of the Riemannian Hessian at the solution, which governs the asymptotic rate of Riemannian gradient descent under the usual local assumptions [15]. Let X⋆ be a critical point of f , and write its eigendecomposition as X⋆ = P diag(d1 , . . . , dn )P ⊤ . With e1 , . . . , en the standard basis of Rn , let ⊤ Eii = P ei e⊤ i P ,
⊤ ⊤ Eij = √12 P (ei e⊤ (i < j) j + ej ei )P
be the frames of Sn built from the eigenbasis of X⋆ , one per unordered pair, so m = n(n + 1)/2 in all. They are orthonormal for the trace inner product. Write B = ∇2 f (X⋆ ) for the Euclidean Hessian at X⋆ , taken positive definite so that condition numbers are defined, and bij = tr Eij B[Eij ] (9) for its diagonal on those frames. Call B a Schur multiplier when it rescales each entry of the eigenbasis without mixing entries, so that B[Eij ] = bij Eij . Proposition 1 (Riemannian Hessian at a critical point). At a critical point X⋆ of f , in every member of (4), Hess f (X⋆ ) = G −1 ◦ B,
(10)
which is self-adjoint for g ϕ and so has real spectrum. Its matrix in the g ϕ -orthonormal frames is B scaled row and column by ϕ(di , dj )1/2 , so its diagonal there is λij = ϕ(di , dj ) bij .
(11)
n ϕ Since Sn ++ is open in S , the Levi-Civita connection of g is
Proof. the directional derivative plus a term bilinear in its two arguments. Differentiating grad f = G −1 [∇f ] therefore splits Hess f (X)[U ] into G −1 ∇2 f (X)[U ] and a remainder of two pieces, the derivative of G −1 applied to ∇f (X) and the connection term evaluated at grad f (X). Both are linear in ∇f (X). At a critical point ∇f (X⋆ ) = 0, so the remainder vanishes and (10) follows. The remainder carried every appearance of the derivative of the metric, which is why only G −1 survives. Self-adjointness follows because g ϕ G −1 [B[U ]], V = tr B[U ]V is symmetric in U and V , since B is a Euclidean Hessian. For the last claim, (7) gives G −1 [Eij ] = ϕ(di , dj )Eij , so the frames stay orthogonal under g ϕ but carry g ϕ (Eij , Eij ) = 1/ϕ(di , dj ), and ϕ(di , dj )1/2 Eij are the g ϕ -orthonormal ones. In these orthonormal frames, (10) gives a rescaled version of the matrix of B in the frames Eij . Each row and column indexed by (i, j) is multiplied by ϕ(di , dj )1/2 . The diagonal entries are therefore ϕ(di , dj )bij , as in (11).
Corollary 2 (Diagonal bound). The condition number κϕ of the Riemannian Hessian Hess f (X⋆ ) obeys κϕ ≥
maxij λij minij λij
(12)
for every B, since a diagonal entry of a symmetric matrix is a Rayleigh quotient. Equality holds when B is a Schur multiplier, where the frames are eigenvectors and the λij are the whole spectrum. The metric reaches the Hessian only through G −1 and never through its derivative, so choosing a member of (4) is always choosing a diagonal preconditioner in the eigenbasis of the current iterate, a valid local choice near the critical point. Write κp,q for κϕ at ϕ = ϕp,q , and κr when that member is the diagonal one p = q = r/2. Plain κ = κ(X⋆ ) = dmax /dmin is the condition number of X⋆ itself. Two particular cases on B are useful to consider: (i) B has the diagonal pure power of degree γ − 2 when bii = c diγ−2 for every i, and (ii) it has the full pure power when all entries of the frame diagonal obey bij = c (di dj )(γ−2)/2 . The second implies the first. Define r⋆ = 2 − γ throughout. Here c > 0 is a scale factor shared by all frames, and it affects no condition number below, since κϕ is a ratio of eigenvalues. Proposition 3 (Conditioning floor). Let ϕ be any weight in (3) that is positively homogeneous of degree r, and let B have the diagonal pure power. Then κϕ ≥
⋆ maxij λij ≥ κ |r−r | . minij λij
(13)
This bound holds independently of the values of bij for i < j and of the couplings between distinct frames. Proof. The first inequality is (12). For the second, only the pairs (i, i) are needed. Homogeneity gives ϕ(d, d) = dr ϕ(1, 1), so λii = ⋆ c ϕ(1, 1) dir+γ−2 , whose ratio at dmax and at dmin is κ r−r . The largest λij over the smallest is therefore at least that ratio, and at least its reciprocal since either of the two may be the larger, hence at least ⋆ κ |r−r | . It should be noted that the floor (right hand side) binds every kernel metric of degree r, not only members of (4). This is true for the Alpha-Procrustes metrics [16]. Their metric operator is diagonal in the frames Eij and equal to d2α−2 on the Eii by [16, Theorem 4], so their weight carries degree 2(1 − α) and (13) puts their floor at ⋆ κ |2(1−α)−r | , which at r⋆ = 0 is the κ 2|α−1| law they report alongside [16, Theorem 6]. The paper [16] recommends α = 1, whose ⋆ degree is zero, so the floor it inherits is κ |r | . For objectives with ⋆ r = 0 this is 1, and (13) leaves the member free. Corollary 4 (When the floor is attained). Assume in addition that B is a Schur multiplier and has the full pure power. Every λij is then the value at (di , dj ) of the single function λ(x, y) = c ϕ(x, y) (xy)(γ−2)/2 ,
x, y > 0.
Furthermore, if λ is monotone in each argument, then both inequalities in (13) are equalities. Proof. The first inequality becomes an equality by Corollary 2. For the second, a symmetric λ monotone in one argument is monotone the same way in the other, so its extremes over the pairs (di , dj ) sit where both arguments are extreme, at (dmax , dmax ) and (dmin , dmin ).
Both extrema occur at diagonal frames Eii . By the proof of Propo⋆ sition 3, the ratio of the maximum to the minimum is κ |r−r | . Thus both inequalities in (13) are equalities. The diagonal member always meets this hypothesis. At p = q = r/2 the function λ is the pure power c (xy)(r+γ−2)/2 , hence monotone at every r. The right-hand side of (13) is minimized at r = r⋆ , where it ⋆ equals one. At other exponents, κϕ is at least κ |r−r | . Minimizing the bound is not the same as minimizing κϕ , and the two coincide when both inequalities become equalities. Whether a given member meets the monotonicity hypothesis is decidable in closed form. Under the full pure power, r⋆ = 2 − γ gives (γ − 2)/2 = −r⋆ /2. Multiplication by this factor shifts both exponents in (4) by −r⋆ /2. The function in Corollary 4 therefore becomes λ(x, y) ∝ xp̃ y q̃ + xq̃ y p̃ ,
⋆
p̃ = p − r2 ,
⋆
q̃ = q − r2 . (14)
The omitted positive factor is common to all frames and cancels from condition numbers. The shifted exponents carry the same two quantities as before, since p̃ + q̃ = r − r⋆ and p̃ − q̃ = p − q. Proposition 5 (Monotonicity criterion). The function (14) is monotone in each argument on (0, ∞)2 if and only if p̃ q̃ ≥ 0, equivalently |p − q| ≤ |r − r⋆ |.
(15)
Proof. The derivative of λ in x is p̃ xp̃−1 y q̃ + q̃ xq̃−1 y p̃ , nonnegative everywhere when p̃, q̃ ≥ 0 and nonpositive everywhere when p̃, q̃ ≤ 0. If instead p̃ > 0 > q̃ then p̃ − 1 > q̃ − 1, so the first term dominates as x → ∞ and the second, which is negative, dominates as x → 0+ , and the derivative changes sign. The case q̃ > 0 > p̃ is the same with the terms exchanged. The second form follows from 4p̃ q̃ = (r − r⋆ )2 − (p − q)2 . How far a member may sit off the diagonal and stay monotone is therefore how far its exponent sits from optimal. The diagonal p = q meets (15) at every exponent, which is why we fix p = q = r/2 and tune r alone. Monotonicity is a sufficient condition for attaining the floor, not a necessary one, and the next proposition specifies by how much. Under the same hypotheses, the conditioning is available in closed form, which settles what happens off the diagonal rather than only when the floor is met. Proposition 6 (Exact conditioning in the pure-power Schur regime). Assume B is a Schur multiplier with the full pure power, and write s = r − r⋆ for the degree error and a = p − q for the shape. Then n o (16) κp,q = max κ |s| , κ |s|/2 cosh a2 log κ . In particular κp,q = cosh( a2 log κ) at r = r⋆ , so once κ > 1 the diagonal member is the unique minimizer and every other member misses the floor by a factor that grows with |p − q|. Proof. Since B is a Schur multiplier, Corollary 2 makes the λij the whole spectrum, so κp,q is the ratio of the largest eigenvalue to the smallest, and Corollary 4 makes each one λ(di , dj ). Up to a positive constant, (14) is (4) with both exponents lowered by r⋆ /2, so the factorization (6) applies to it with s in place of r and a in place of p − q, a x λ(x, y) ∝ (xy)s/2 cosh log , (17) 2 y
the separation the family was built on, now read on the spectrum. What remains is the largest and the smallest of (17) over the frames, and a logarithm puts both within reach, since it turns that product into a sum. Write ℓi = log di . At a frame Eij the logarithm of (17) is an affine function of (ℓi , ℓj ) carrying s alone, plus log cosh(a(ℓi − ℓj )/2), which is convex and nonnegative. The values below are scaled by (dmin dmax )−s/2 , which the ratio cancels. p Minimum. The affine part is least over the frames where di dj is smallest if s ≥ 0 and largest if s < 0, and either end forces di = dj , where the convex part vanishes. The two are least together, at κ−|s|/2 . Maximum. The sum is convex on the square [ℓmin , ℓmax ]2 , so it is maximal at a corner, and every corner is a frame. The two diagonal corners are the Eii at dmin and at dmax , with values κ−s/2 and κ s/2 , and the other two are the single frame pairing dmin with dmax , where the power in (17) cancels that scaling exactly and the value is cosh( a2 log κ). The maximum is therefore max κ |s|/2 , cosh( a2 log κ) , and dividing it by the minimum gives (16). Only κ enters, not the interior of the spectrum. By (16) a member meets the floor κ |s| exactly when cosh( a2 log κ) ≤ κ |s|/2 , that is when log 2 |p − q| ≤ log2 κ arccosh κ |s|/2 = |r − r⋆ | + 2log + O κ−|s| . κ (18) This is wider than the monotonicity budget (15) by 2 log 2/ log κ, so Proposition 5 is the limit of (18) as κ grows, and the slack it leaves out is the room a non-monotone member has to attain the floor anyway. Inside that budget (16) returns κ |s| whatever the shape is, so κp,q is flat in p − q across a band and grows like κ |p−q|/2 only outside it. The band closes exactly at r = r⋆ , where (18) has width zero. Within the pure-power Schur regime, the diagonal member p = q = r/2 minimizes κp,q at each fixed r. Thus it suffices to tune r. The same reading weighs the two parameters against each other: moving the degree by t costs κ t , while moving the shape by t costs cosh( 2t log κ), which carries half that exponent, so the degree is worth twice the shape. 4. SELECTION RULES FOR r, p AND q In the pure-power Schur regime, Corollary 4 settles p and q once r is fixed, namely p = q = r/2. That member attains (13) at every r, and at r = r⋆ it is the only member (15) admits. It is also the cheapest member, by (5), and it names the exponent as r⋆ = 2 − γ, so only γ remains to be estimated. A Hessian that is a power congruence. Proposition 3 constrains the frame diagonal, so a rule must be stated in those terms. The condition is that the Hessian act by congruence with a power of X, ∇2 f (X)[U ] = c X (γ−2)/2 U X (γ−2)/2
(19)
for some c > 0, at X = X⋆ , since that gives bij = c (di dj )(γ−2)/2 exactly, which is the full pure power, and r⋆ = 2 − γ. Fix T and C in n S++ . The Hessian of 12 kX − T k2F is U 7→ U , so γ = 2 and r⋆ = 0. The Hessian of tr(CX) − log det X is U 7→ X −1 U X −1 , so γ = 0 and r⋆ = 2. The Hessian of 21 kX −1 − T k2F is U 7→ X −2 U X −2 at its minimizer X⋆ = T −1 . Thus γ = −2 and r⋆ = 4, above the degrees 0, 1 and 2 of the Euclidean, Bures–Wasserstein and affineinvariant metrics, respectively. The first two meet (19) at every X,
so they name r⋆ before X⋆ is known, while the third meets it only at the minimizer. A Hessian of mixed degree. A sum of terms of different degrees, such as an evidence lower bound or a regularized loss, meets no single (19), so thepexponent has to be estimated rather than read off. Write uij = log di dj . The diagonal member has ϕ(di , dj ) = (di dj )r/2 = e ruij by (5), so (11) reads log λij = r uij + log bij , affine in r. The spread of the frame diagonal in the logarithm, which is the spectrum itself only when B is a Schur multiplier, is then a quadratic in r with a closed-form minimizer, P X 2 i≤j (uij − ū) log bij P , r̂ = arg min log λij − log λ = − 2 r i≤j (uij − ū) i≤j
(20) where a bar is the mean over the m frames. The right-hand side is the negative of the least-squares slope of log bij against uij . Thus r̂ can be computed in a single pass over the frame data. The estimate is well defined when all bij > 0 and the uij are not all equal. The latter condition fails precisely when X = cI for some c > 0. In the pure power case, it recovers r⋆ exactly. From here we index the m frames by k = 1, . . . , m and write uk and vk = log bk for the pair each one carries, so that log λk = r uk + vk . The sum of squares is a surrogate. What (12) depends on is the range of log λk rather than its spread, and minimizing that range, ř = arg min R(r),
r R(r) = max r uk + vk − min r uk + vk ,
(21)
k≤m
k≤m
is a one-dimensional convex problem on the same m numbers, since a maximum of affine functions is convex. Both criteria are minus the slope of v against u, (20) fitted in ℓ2 and (21) in ℓ∞ . The logarithm of (12) is a range, and R(r) = 2 minc maxk |r uk + vk − c| depends on the two extreme frames alone, where the sum of squares depends on all m. Note this can be solved as a linear program [17]. A sampling approach to compute ř and r̂ efficiently. Both criteria read the same m numbers bij , one per frame, a count quadratic in n. A single eigendecomposition of X supplies the di and the frames, and each bij of (9) then costs one Hessian-vector product, since it pairs Eij against ∇2 f (X)[Eij ], which a directional derivative of the gradient delivers without ever forming B. Tuning is therefore one eigendecomposition and m products, paid once against a run of many iterations. A slope is a two-point quantity, though, and the frames are far from equally informative about it, so it is worth asking what a subset costs. Proposition 7 (Estimation with fewer frames). Fix X ∈ Sn ++ with κ = dmax /dmin > 1 and let uk , vk for k = 1, . . . , m be its frame data, so that log λk = r uk + vk on the diagonal member by (11) and (5). Write S0 for the pair of diagonal frames Eii at di = dmin and at di = dmax . Since uij = 12 (uii + ujj ), every uk lies in [log dmin , log dmax ], and S0 attains both endpoints. Let r̂ and ř minimize (20) and (21) over all m frames, and let ε ∈ Rm be the residual of the fit, εk = (vk − v̄) + r̂ (uk − ū). (22) For a nonempty S ⊆ {1, . . . , m} write r̂S and řS for the minimizers of (20) and (21) taken over S alone, let εS P and RS be the residual and the range restricted to S, and put σS2 = k∈S (uk − ūS )2 with ūS the mean of u on S. Then the following hold. (i) If σS > 0, then r̂S − r̂ ≤
kεS k2 . σS
(23)
(ii) If S ⊇ S0 , then 0 ≤ R(řS ) − R(ř) ≤ 4kεk∞ .
(24)
Proof. (i) Restricted least squares on S is r̂S = −
1 X (uk − ūS ) vk . σS2
objective
k∈S
Substitute vk = v̄ − r̂(uk − ū) P+ εk , which is (22) rearranged. The two constants vanish against k∈S (uk − ūS ) = 0 and the linear term returns r̂, so r̂S = r̂ −
1 X (uk − ūS ) εk , σS2
and Cauchy–Schwarz bounds that sum by σS kεS k2 . (ii) By (22), log λk = (r − r̂)(uk − ū) + εk up to a constant, which no range sees. A range of a sum is at most the sum of the ranges, so R(r) ≤ |r − r̂| log κ+2kεk∞ , while RS (r) ≥ RS0 (r) ≥ |r − r̂| log κ − 2kεk∞ because u spans log κ already on S0 and enlarging a set only widens its range. Hence R − RS ≤ 4kεk∞ at every r. Dropping frames lowers a maximum and raises a minimum, so RS ≤ R pointwise and R(ř) ≥ RS (ř) ≥ RS (řS ). Subtracting that from R(řS ) bounds it by R(řS ) − RS (řS ) and so by 4kεk∞ , and ř minimizes R, which gives the first inequality.
r̂S0 − r̂ log κ ≤ 2kεk∞ .
r⋆ Eucl.
BW aff.-inv. log-Eucl.
least squares 0 0 3.99 precision 2 7.99 3.99 inverse fit 4 15.97 11.98
7.99 0 7.99
7.99 2.07 7.99
AP
tuned
0.30 7.99 15.97
0 0 0
5.1. The exponent read off the Hessian
k∈S
Note that on the extreme pair alone σS2 0 = √ kεS0 k2 ≤ 2 kεk∞ , so (23) leads to
Table 1. The congruence rule on a real covariance, n = 36 and κ(X⋆ ) = 9836. Entries are log10 κϕ , so (13) predicts 3.99 |r − r⋆ | from the degree r of the weight alone. The three objectives are those of Section 4 in order, and AP is the Alpha-Procrustes member at α = 1. The tuned members in the first two rows are Euclidean (r⋆ = 0) and affine-invariant (r⋆ = 2), respectively.
1 log2 κ 2
and
(25)
An objective built from a distance. For a squared-distance objective, we use the metric that defines the distance. The Riemannian Hessian of 12 dist2 (·, T ) at X⋆ = T is the identity in the metric that defines dist, so that metric is exactly optimal for it. We call this reciprocity. A least-squares fit 21 kX − T k2F corresponds to the Euclidean metric, and 12 dist2BW (·, T ) to Bures–Wasserstein, which the fitted exponent cannot reach because it sits off the diagonal p = q. Likewise 21 dist2AI (·, T ) corresponds to affine-invariant geometry, and 1 klog X − log T k2F to log-Euclidean, a weight the family does not 2 contain at all. A Fréchet mean averages several such terms. Under the Euclidean and log-Euclidean metrics the Hessian is still the identity at the barycenter, so the match stays exact, while under the others the minimizer is none of the Ti and the match is only as close as the spread of the data allows. This observation is consistent with (19). By (5), that condition makes the Euclidean Hessian a positive multiple of the metric operator at p = q = (2 − γ)/2. It therefore gives the same local Hessian matching within the diagonal family. 5. EXPERIMENTS The experiments below score the tuned member against the named metrics a practitioner would otherwise pick. Every run below uses Algorithm 1 with the retraction (8) and an Armijo line search [18, 19], starts from the same X0 and stops at the same tolerance, so only the weight differs. The tolerance is a relative objective gap. Conditioning is always κϕ at a common reference solution.
We minimize the three objectives of Section 4, the least-squares fit, the Gaussian precision estimate and the inverse fit, whose powercongruence Hessians name r⋆ = 0, 2 and 4 before anything is run. Here T = C is a real covariance, the class-conditional pixel covariance of a standard digits set at n = 36, and all three problems share κ(X⋆ ) = 9836 because κ(A) = κ(A−1 ), so the objective alone moves r⋆ . A power congruence is a full pure power and a Schur multiplier, so Proposition 6 gives κp,q in closed form, and for ⋆ a diagonal weight of degree r it reduces to κ |r−r | , which Table 1 confirms to the digit. Off the diagonal a weight of the right degree need not attain the floor, and the table shows both outcomes. Bures– Wasserstein does attain it, because the extreme λij fall on the frames Eii where its weight is a multiple of the diagonal one of the same degree. It does so in all three rows only because no r⋆ among them lies near its own degree. Read at (p, q) = (1, 0), (16) meets the floor exactly when |1 − r⋆ | ≥ log2 κ log cosh( 12 log κ), which is 0.85 here and rises toward 1 as κ grows, against the 1, 1 and 3 the three rows supply. At r⋆ = 1 its degree is exactly right and it still pays cosh( 12 log κ) = 49.6√where the diagonal member p = q = 21 pays 1, a penalty of order κ/2. Matching the degree is therefore never enough on its own. Log-Euclidean is the other outcome, carrying the right degree on the second objective and still missing the floor by 116 because it is built from the logarithmic rather than the geometric mean, so the degree is necessary and not sufficient. The third row is the one to note. A standard estimator has r⋆ = 4, two degrees past affine-invariant, and there the best named weight is 9.7 × 107 worse conditioned than the diagonal member at r⋆ . The frame diagonal is a pure power here, so vk is exactly affine in uk , the residual of Proposition 7 vanishes, and its pair S0 suffices: two Hessian-vector products in place of m = 666 return each of r⋆ = 0, 2 and 4 to twelve digits. 5.2. An exponent past every named metric The objective is a Mahalanobis matrix learned from labeled triplets under a regularizer [20, 21], f (X) =
1 X s 1 + hX, Dt i + µR(X), |T | t∈T
(26)
where s(z) = log(1 + ez ) is the softplus and Dt = (xa − x+ )(xa − x+ )⊤ − (xa − x− )(xa − x− )⊤ for a triplet of an anchor xa , a positive x+ and a negative x− . We use the Wine and Breast Cancer Wisconsin datasets from the UCI repository [22], as distributed with scikit-learn [23], with a stratified 70/30 split repeated three times, |T | = 2000 triplets per split, µ = 0.01 and X0 = I. Features are centered. We consider three regularizers: R = tr X − log det X,
tr(X 1)
101 10 4 10 9
1 tr(X 2) 2
/ best
Wine f f
10 4 10 9
Breast f f
logdet X
range criterion
least-squares fit
101
101
2 1 102
frames read
0
50
Eucl. BW affine-inv.
100 0
50
iterations
log-Eucl. AP, = 1
100 0
r, L = 5 r, L = 5
50
100
r, L = 25 r, L = 25
Fig. 1. Objective gap for (26) against iterations at µ = 0.01. Rows correspond to the two datasets, and columns to the three regularizers. The regularizers are labeled by R − tr X, as in Table 2. All runs start from X0 = I and use the same line search. The two-stage curves run at the default r = 2, the affine-invariant weight, until the marker at L and branch there. A five-step pilot suffices on the barriers but not under the log-determinant, which is why L = 25 is used throughout. Log-Euclidean carries open squares because it can sit under affineinvariant to within half a percent.
tr X + tr(X −1 ), and tr X + 12 tr(X −2 ). Their diagonal frame cur−3 vatures are d−2 and 3d−4 i , 2di i , respectively. For the regularizers alone, these correspond to r⋆ = 2, 3 and 4. We call the last two barriers. The log-determinant is the usual choice on Sn ++ , and the exponent it sets is exactly affine-invariant, so nothing is left to tune there. The two barriers grow faster at the boundary and carry r⋆ past every named metric, which is what makes them worth running. The linear term adds nothing to the Hessian and is there only to bound X⋆ above. The triplet loss carries no degree at all, so the sum carries none, no rule of Section 4 names the exponent and it has to be estimated, from the least-squares fit (20) or the range criterion (21). Both read the frame diagonal bij and the eigenvalues di at a point, and the point that matters is the solution being sought, so we run Algorithm 1 in two stages. Descend L steps from X0 at the default exponent r = 2, evaluate the criterion at the point reached, then continue from that point at the exponent it returns. Every iteration count we report includes the L pilot steps. Reading the frame diagonal costs m Hessian-vector products in general, but they batch into one pass over the triplets here, so either criterion costs 0.9 to 1.5 iterations. Figure 1 plots the objective gap against iterations for all six settings. On the four barrier panels both criteria reach machine precision inside the budget and no named weight does, because the exponent they return sits near 3 or 4 and the nearest named degree is 2. The range criterion returns 3.45 and 4.15 on Wine and 3.10 and 4.05 on Breast under the two barriers, where the least-squares fit runs higher, to 4.87. Under the log-determinant both return near 2, and that column is where the criteria add least, since affine-invariant already carries that regularizer’s exponent. Figure 2 varies the number of frames used for fitting from two to m. Each subset contains the pair S0 from Proposition 7. The remaining frames are sampled at random. Dropping frames costs little, and for one of the two criteria it gains. The range criterion is flat, so the extra frames are not required. The least-squares fit, on the
Wine, logdet X Wine, tr(X 1) Wine, 12 tr(X 2)
102
frames read
Breast, logdet X Breast, tr(X 1) Breast, 12 tr(X 2)
Fig. 2. Conditioning delivered against the number of frames read, over the six settings of (26) at µ = 0.01, with the range criterion (21) on the left and the least-squares fit (20) on the right. Curves are named by dataset and by R − tr X as in Table 2. Each value is κϕ at the exponent the criterion returns, divided by the best any diagonal member attains, so 1 is the floor of the sweep. Each subset contains the pair S0 from Proposition 7. Additional frames are sampled at random. The curves show the median over twelve draws at each subset size. Reading more frames does nothing for the range criterion and hurts the least-squares fit.
other hand, is better off with fewer. No subset drawn violated (23), (24) or (25). Read with Table 1, where the residual vanishes so that any two Hessian-vector products reproduce every r⋆ exactly, the n diagonal frames suffice in both regimes. Table 2 sweeps the regularization weight, since µ = 0.01 is a single setting and the advantage need not survive elsewhere. Each entry is a gain, the condition number of the best named metric divided by the one the criterion delivers, so a value above one is the factor by which tuning improves the conditioning and a value below one is a loss. The two barriers gain in every setting, and by the widest margin where the regularization is weakest. The log-determinant is the exception, gaining little and sometimes losing, since affine-invariant already carries the exponent it sets. On Alpha-Procrustes. The weights of [16] carry degree 2(1 − α), so the member at α = 1 they recommend has degree zero, like the Euclidean case. Its weights are given by 2(x2 + y 2 )/(x + y)2 . This equals 1 at x = y and approaches 2 as x/y tends to 0 or ∞. It agrees with the Euclidean weight on the diagonal frames Eii , but differs on off-diagonal frames with di = 6 dj . Table 1 compares the resulting Hessian condition numbers. For the precision and inversefit objectives, the extreme Hessian eigenvalues occur on the diagonal frames, so the two metrics have the same condition number. For the least-squares objective, B is the identity. Euclidean then has condition number 1, while the variation of the Alpha-Procrustes weight gives a condition number (close to) 2. Figure 1 also shows similar convergence curves for these two metrics from X0 = I. 5.3. Tuning the shape We learn an SPD task covariance for the seven torque outputs in the SARCOS inverse-dynamics dataset [24]. A multi-output Gaussian process uses this matrix to model dependence between the outputs [25]. Let Y ∈ RN ×7 contain the standardized torques and let y = vec(Y ). For an input covariance matrix Γ ∈ RN ×N and noise variance ν, the model has Σ(X) = X ⊗ Γ + νI7N ,
X ∈ S7++ .
(27)
ř
r̂
µ = 0.01
ř
ř
r̂
r̂
Wine − log det X 1.06 1.02 1.67 1.19 tr(X −1 ) 1.68 1.67 4.20 3.90 −2 1 tr(X ) 1.96 1.92 5.54 4.68 2
1.00 0.74 10.8 10.2 17.1 11.4
Breast − log det X 1.04 1.02 1.05 1.04 tr(X −1 ) 2.63 2.36 6.43 4.38 −2 1 tr(X ) 3.88 3.29 10.9 6.55 2
0.92 0.98 14.3 7.32 31.1 11.8
q p=
δ ⋆ p=
δ ⋆
q−
q= p−
p
0
Fit r0
Choose δ ⋆
Refit r1
Fig. 3. Selecting the exponent and shape in the (p, q) family. Blue marks the first exponent fit on p = q. Yellow selects δ ⋆ along p+q = r0 , and green refits the exponent along p − q = δ ⋆ . The gray branch contains the same metrics with p and q exchanged. The positions are schematic. The arrows give the selection order. The final fit can move in either direction along the green line.
We minimize the Gaussian negative log marginal likelihood [26], i 1 h f (X) = log det Σ(X) + y⊤ Σ(X)−1 y . (28) 2N The input covariance Γ uses a radial basis function (RBF) kernel. We calibrate its parameters and ν once at X = I, then hold them fixed across methods. We use seven training subsets, four with N = 512, two with N = 2,048 and one with N = 4,096. We standardize each training subset separately. As in Section 5.2, we descend L = 25 steps from X0 = I at the affine-invariant exponent r = 2, then tune at the point reached XL . The preceding experiments fit r on the diagonal p = q. Here we also tune the shape a = p−q of Proposition 6. Exchanging p and q leaves the metric unchanged, so we work with its magnitude δ = |a|. The separation of exponent and shape suggests a tuning heuristic based on the frame data at XL . After the first exponent fit, we use the shape to raise the smaller entries of the scaled diagonal toward their maximum. We then refit the exponent using the adjusted diagonal. At XL , read the frame diagonal bij and the eigenvalues di as in Section 4. All seven pilot points have positive frame curvatures and pass the numerical non-degeneracy checks for exponent fitting. Keep the frame data uk , vk , with bk = bij for the frame associated with (i, j). Write ξk = 21 log(di /dj ) for that pair and evaluate the scaled diagonal (11) at XL . The factorization (6) gives, after removing the 2 common factor 4δ , ek (r, δ) = 4−δ2 λk (r, δ) = eruk +vk cosh(δξk ). λ
r0
R − tr X
µ = 0.1
q=
µ=1
q
p+
Table 2. Conditioning gains from tuning the exponent for differ ent regularization weights. Each entry is minϕ κϕ /κr , where the minimum is taken over the Euclidean, Bures–Wasserstein and affineinvariant metrics. We select r after the L = 25 pilot using either the range criterion ř in (21) or the least-squares fit r̂ in (20). Both criteria use all m frames. A value above one gives the factor by which tuning improves the conditioning. The two criteria agree closely at µ = 1 and separate as the regularizer weakens, with ř generally ahead.
(29)
The factor cancels from the ratio of the largest to the smallest entry. At δ = 0, let r0 be r̂ from the least-squares fit (20) or ř from the ek (r0 , 0). Choose the largest range criterion (21). Set W0 = maxk λ ek (r0 , δ) is at most W0 . For a nonscalar δ ≥ 0 for which every λ spectrum, this gives ek (r0 , 0) arccosh W0 /λ δ ⋆ = min . (30) ξk >0 ξk Each frame with ξk > 0 gives an upper bound on δ. For 0 ≤ δ ≤ δ ⋆ , the maximum normalized entry stays at W0 and the minimum is nondecreasing. Thus (30) minimizes their ratio over this interval.
To obtain r1 , repeat the fitting rule from the first step with δ = δ ⋆ in (29). The selected metric has p=
r1 + δ ⋆ , 2
q=
r1 − δ ⋆ . 2
(31)
We select these parameters once and continue descent from XL . In Figure 4, Shape-LS uses the least-squares fit and Shape-Range uses the range criterion. Figure 3 shows these three steps in the parameter plane. Geometrically, this three-step construction can reach any member of the (p, q) family, up to exchanging p and q. We use the diagonal members returned by the same fitting rules as the r-only baselines. We also compare with Euclidean, BW, affineinvariant and log-Euclidean. All eight methods continue from XL with the same Armijo history. The four fitted methods read all m = 7(7 + 1)/2 = 28 entries bij of the frame diagonal. Runs stop when [f (X) − f (Xref )]/[f (I) − f (Xref )] ≤ 10−8 , with a total budget of 2,000 iterations. Here Xref is a common numerical reference solution. We evaluate κϕ from the full preconditioned Hessian at Xref . With either fitting rule, the selected shape is positive in six cases (Fig. 4(a)). The remaining case, N = 512 with seed 3, gives δ ⋆ = 0. The full condition number improves over the corresponding diagonal baseline in the six positive-shape cases and stays unchanged in the seventh. Median κϕ falls from 7.70 to 3.53 for the least-squares fit and from 7.64 to 3.28 for the range criterion (Fig. 4(c)). Relative to the corresponding diagonal member, the median reduction is 55.80% for the least-squares fit and 54.54% for the range criterion (Fig. 4(d)). Both choices need a median of 28 iterations, including the pilot. Each diagonal member needs 31, log-Euclidean needs 30 and affine-invariant needs 33. Both fitting rules read the diagonal bij of the Euclidean Hessian in the frames Eij . At a critical point where B is a Schur multiplier, the entries λij give the whole Riemannian Hessian spectrum by Corollary 2. Coupling between frames makes the diagonal ratio a surrogate for κϕ . The diagonal ratio is non-increasing under the range refit. The least-squares refit minimizes the variance of the log diagonal. Shape tuning reuses di , the frames and the diagonal bij from the first exponent fit. It adds one exponent fit and an O(m) pass
Eucl. BW
affine-inv. log-Eucl.
r-only LS
(a) Selected parameters
(b) Convergence after the pilot
LS Range (r0, 0) (r1, δ ⋆ )
0.4 0.2 0.0 r-only subfamily: δ = 0
1.70
1.75
1.80
10−5
Gap (×10−8)
0.6
Zoom near iteration 2
Relative objective gap
Shape δ = |p − q|
0.8
Shape-LS Shape-Range
r-only Range
10−6
10−8 0
1
102 10 2 affine- log- r-only r-only Shape- Shapeinv. Eucl. LS Range LS Range
Condition number reduction (%)
Condition number κϕ
103
BW
2
3
4
5
6
7
8
9
Iterations after the pilot
One case Median
Eucl.
3.8
10−7
1.85
(c) Hessian conditioning
10
4.0
1.95 2.00 2.05 Iterations after pilot
Exponent r = p + q
4
4.2
(d) Gain over the r-only baseline Median: 55.80% (LS), 54.54% (Range)
60
40
20
0 512 s0
512 s1
512 s2
512 s3
2,048 2,048 4,096 s1 s2 s3
Training size / seed
Fig. 4. Exponent and shape tuning on SARCOS over seven cases. (a) Open markers show (r0 , 0) and filled markers show (r1 , δ ⋆ ). (b) Median relative objective gap after the shared 25-iteration affine-invariant pilot. Completed runs are held at 10−8 for aggregation. The inset enlarges the second iteration after the pilot. (c) Condition number κϕ of the full preconditioned Hessian at Xref . Each point is a case and black bars mark medians. (d) Percentage reduction 100(1 − κϕ /κr0 ) relative to the corresponding r-only baseline. Each baseline uses the initial exponent r0 from the same fitting rule with δ = 0.
for (30). No additional Hessian-vector products are needed. Both least-squares fits are closed form, so the arithmetic after reading bij remains O(m). For the range criterion, we solve two linear programs instead of one. The fitting data still require O(m) storage. 6. CONCLUSION The family of metrics defined by X p LX q + X q LX p = U contains Euclidean, Bures–Wasserstein and affine-invariant geometry as exact members and reaches past all three. Each iteration of Algorithm 1 uses an eigendecomposition and the retraction (8), which preserves positive definiteness. The conditioning of the Riemannian Hessian at the solution obeys a floor set by r = p + q alone, and in the purepower Schur regime the diagonal member p = q = r/2 attains that floor, so the two-parameter choice collapses to one number. We give principled ways to estimate r, and show where tuning it helps and where a named metric already suffices. The SARCOS example also shows a further gain from tuning the shape.
7. REFERENCES [1] Xavier Pennec, Pierre Fillard, and Nicholas Ayache, “A Riemannian framework for tensor computing,” International Journal of Computer Vision, vol. 66, no. 1, pp. 41–66, 2006. [2] Alexandre Barachant, Stéphane Bonnet, Marco Congedo, and Christian Jutten, “Multiclass brain–computer interface classification by Riemannian geometry,” IEEE Transactions on Biomedical Engineering, vol. 59, no. 4, pp. 920–928, 2012. [3] Asuka Takatsu, “Wasserstein geometry of Gaussian measures,” Osaka Journal of Mathematics, vol. 48, no. 4, pp. 1005–1026, 2011. [4] Luigi Malagò, Luigi Montrucchio, and Giovanni Pistone, “Wasserstein Riemannian geometry of Gaussian densities,” Information Geometry, vol. 1, no. 2, pp. 137–179, 2018. [5] Rajendra Bhatia, Tanvi Jain, and Yongdo Lim, “On the Bures–
Wasserstein distance between positive definite matrices,” Expositiones Mathematicae, vol. 37, no. 2, pp. 165–191, 2019. [6] Vincent Arsigny, Pierre Fillard, Xavier Pennec, and Nicholas Ayache, “Geometric means in a novel vector space structure on symmetric positive-definite matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 29, no. 1, pp. 328–347, 2007. [7] Yann Thanwerdas and Xavier Pennec, “O(n)-invariant Riemannian metrics on SPD matrices,” Linear Algebra and its Applications, vol. 661, pp. 163–201, 2023. [8] Fumio Hiai and Dénes Petz, “Riemannian metrics on positive definite matrices related to means,” Linear Algebra and its Applications, vol. 430, no. 11-12, pp. 3105–3130, 2009. [9] Yann Thanwerdas and Xavier Pennec, “Is affine invariance well defined on SPD matrices? a principled continuum of metrics,” in Geometric Science of Information. 2019, pp. 502–510, Springer. [10] Yann Thanwerdas and Xavier Pennec, “The geometry of mixed-Euclidean metrics on symmetric positive definite matrices,” Differential Geometry and its Applications, vol. 81, pp. 101867, 2022. [11] Bamdev Mishra and Rodolphe Sepulchre, “Riemannian preconditioning,” SIAM Journal on Optimization, vol. 26, no. 1, pp. 635–660, 2016. [12] Andi Han, Bamdev Mishra, Pratik Jawanpuria, and Junbin Gao, “On Riemannian optimization over positive definite matrices with the Bures–Wasserstein geometry,” in Advances in Neural Information Processing Systems, 2021, vol. 34. [13] Andi Han, Bamdev Mishra, Pratik Jawanpuria, and Junbin Gao, “Learning with symmetric positive definite matrices via generalized Bures–Wasserstein geometry,” in Geometric Science of Information. 2023, pp. 405–415, Springer. [14] Boris Shustin and Haim Avron, “Riemannian optimization with a preconditioning scheme on the generalized Stiefel manifold,” Journal of Computational and Applied Mathematics, vol. 423, pp. 114953, 2023. [15] Bin Gao, Renfeng Peng, and Ya-xiang Yuan, “Optimization on product manifolds under a preconditioned metric,” SIAM Journal on Matrix Analysis and Applications, vol. 46, no. 3, pp. 1816–1845, 2025. [16] Derun Zhou, Keisuke Yano, and Mahito Sugiyama, “Riemannian optimization over symmetric positive definite matrices with the Alpha-Procrustes geometry,” 2026, arXiv:2605.00396. [17] Stephen Boyd and Lieven Vandenberghe, Convex Optimization, Cambridge University Press, 2004. [18] P.-A. Absil, Robert Mahony, and Rodolphe Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, Princeton, NJ, 2008. [19] Nicolas Boumal, An Introduction to Optimization on Smooth Manifolds, Cambridge University Press, 2023. [20] Kilian Q. Weinberger and Lawrence K. Saul, “Distance metric learning for large margin nearest neighbor classification,” Journal of Machine Learning Research, vol. 10, pp. 207–244, 2009. [21] Jason V. Davis, Brian Kulis, Prateek Jain, Suvrit Sra, and Inderjit S. Dhillon, “Information-theoretic metric learning,” in Proceedings of the 24th International Conference on Machine Learning, 2007, pp. 209–216.
[22] Dheeru Dua and Casey Graff, “UCI machine learning repository,” University of California, Irvine, School of Information and Computer Sciences, 2019. [23] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay, “Scikitlearn: Machine learning in Python,” Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011. [24] Sethu Vijayakumar and Stefan Schaal, “Locally weighted projection regression: Incremental real time learning in high dimensional space,” in ICML ’00 Proceedings of the Seventeenth International Conference on Machine Learning. 2000, pp. 1079–1086, Morgan Kaufmann Publishers Inc. [25] Edwin Bonilla, Kian Chai, and Christopher Williams, “Multitask Gaussian process prediction,” in Advances in Neural Information Processing Systems, J. Platt, D. Koller, Y. Singer, and S. Roweis, Eds. 2007, vol. 20, Curran Associates, Inc. [26] Carl Edward Rasmussen and Christopher K. I. Williams, Gaussian Processes for Machine Learning, MIT Press, Cambridge, MA, 2006.