Semidefinite approximations of the matrix logarithm

TL;DR

This paper develops semidefinite approximation methods for the matrix logarithm that preserve its operator concavity, enabling efficient convex optimization involving this function in quantum information and related fields.

Contribution

It introduces a novel approach to approximate the matrix logarithm with semidefinite programs, facilitating practical optimization in quantum information theory.

Findings

Allows use of standard SDP solvers for log-based convex problems

Provides small-size semidefinite representations for operator concave functions

Extends to other operator concave and monotone functions

Abstract

The matrix logarithm, when applied to Hermitian positive definite matrices, is concave with respect to the positive semidefinite order. This operator concavity property leads to numerous concavity and convexity results for other matrix functions, many of which are of importance in quantum information theory. In this paper we show how to approximate the matrix logarithm with functions that preserve operator concavity and can be described using the feasible regions of semidefinite optimization problems of fairly small size. Such approximations allow us to use off-the-shelf semidefinite optimization solvers for convex optimization problems involving the matrix logarithm and related functions, such as the quantum relative entropy. The basic ingredients of our approach apply, beyond the matrix logarithm, to functions that are operator concave and operator monotone. As such, we introduce strategies for constructing semidefinite approximations that we expect will be useful, more generally, for studying the approximation power of functions with small semidefinite representations.

Figures

Tables

Table 1. Table 1: List of functions available in the package CvxQuad . The last column gives the size of the semidefinite representations (here LMI stands for Linear Matrix Inequality, and corresponds to a constraint of the form in ( 1 )). The parameters m 𝑚 m and k 𝑘 k control the accuracy of the approximation (see Proposition 1 ) and n 𝑛 n is the size of the matrix arguments.

op_rel_entr_epi_cone	$-X^{1/2}\log \left(X^{-1/2}Y{X}^{-1/2}\right)X^{1/2}⪯T$	$m+k$ LMIs of size $2n\times 2n$ each
quantum_entr	$\rho \mapsto -\mathrm{Tr}\left[\rho \log \rho \right]$ (Concave)	$m+k$ LMIs of size $2n\times 2n$ each
trace_logm	$\rho \mapsto \mathrm{Tr}\left[\sigma \log \rho \right]$ ($\sigma ⪰0$ fixed; Concave)	$m+k$ LMIs of size $2n\times 2n$ each
quantum_rel_entr	$(\rho ,\sigma )\mapsto \mathrm{Tr}\left[\rho \left(\log \rho -\log \sigma \right)\right]$ (Convex)	$m$ LMIs of size $\left(n^2+1\right)\times \left(n^2+1\right)$
		and $k$ LMIs of size $2n^2\times 2n^2$ each

Table 2. Table 2: Maximum entropy optimization ( 33 ) via our method and the successive approximation scheme of CVX on different random instances. We see that our method can be much faster than the successive approximation method while having the same accuracy. The accuracy of the different methods is measured via the difference | p − p Mosek | 𝑝 subscript 𝑝 Mosek |p-p_{\textsc{Mosek}{}}| where p Mosek subscript 𝑝 Mosek p_{\textsc{Mosek}{}} is the optimal value returned by the built-in Mosek solver for maximum entropy problems ( mskenopt ), and p 𝑝 p is the optimal value returned by the considered approximation method. For our method we used the parameters ( m , k ) = ( 3 , 3 ) 𝑚 𝑘 3 3 (m,k)=(3,3) . The last column also gives the gap between the optimal value returned by the two different approximation methods.

		Successive approximation (CVX)		Padé approximation (this paper)
$n$	$\ell$	time (s)	accuracy	time (s)	accuracy	$\left|p_{sa}-p_{Pade}\right|$
50	25	0.34 s	1.065e-05	0.32 s	1.719e-06	8.934e-06
100	50	0.52 s	1.398e-06	0.34 s	2.621e-06	1.222e-06
200	100	1.10 s	6.635e-06	0.88 s	2.767e-06	3.868e-06
400	200	3.38 s	2.662e-05	0.72 s	1.164e-05	1.498e-05
600	300	9.14 s	2.927e-05	1.84 s	2.743e-05	1.843e-06
1000	500	52.40 s	1.067e-05	3.91 s	1.469e-04	1.362e-04

Table 3. Table 3: Geometric programming ( 34 ) using our method (with ( m , k ) = ( 3 , 3 ) 𝑚 𝑘 3 3 (m,k)=(3,3) ) and the successive approximation scheme of CVX, on different random instances. The column “ sp ” indicates the sparsity of the power vectors a j , k subscript 𝑎 𝑗 𝑘 a_{j,k} (i.e., how many variables appear in each monomial terms). Also we used w 0 = w 1 = ⋯ = w ℓ = 5 subscript 𝑤 0 subscript 𝑤 1 ⋯ subscript 𝑤 ℓ 5 w_{0}=w_{1}=\dots=w_{\ell}=5 (i.e., the posynomial objective as well as the posynomial constraints all have 5 terms). Accuracy is measured via absolute error between the optimal value returned by the approximation and the built-in Mosek solver for geometric programs ( mskgpopt ).

			Successive approximation (CVX)		Padé approximation (this paper)
$n$	$\ell$	sp	time (s)	accuracy	time (s)	accuracy	$\left|p_{sa}-p_{Pade}\right|$
50	50	0.3	1.28 s	2.509e-07	0.94 s	2.106e-06	1.856e-06
50	100	0.3	1.78 s	2.045e-05	1.03 s	3.122e-05	1.077e-05
100	100	0.1	1.57 s	4.759e-06	1.16 s	5.197e-06	4.383e-07
100	150	0.1	3.60 s	8.484e-06	1.60 s	2.240e-06	6.244e-06
100	200	0.1	7.60 s	1.853e-06	2.69 s	3.769e-06	1.916e-06
200	200	0.1	7.47 s	2.441e-07	3.72 s	7.505e-07	9.945e-07
200	400	0.1	42.71 s	3.666e-06	14.36 s	2.855e-06	6.521e-06
200	600	0.1	184.33 s	7.899e-06	35.45 s	4.480e-06	3.419e-06

Table 4. Table 4: Result of solving the optimization problem ( 35 ) for different Hermitian positive definite matrices Y 𝑌 Y of size n × n 𝑛 𝑛 n\times n with Tr ⁡ [ Y ] = 1 Tr 𝑌 1 \operatorname{Tr}[Y]=1 . The problems were implemented using CVX as shown above and solved using SDPT3. The accuracy column reports the quantity | p − 1 | 𝑝 1 |p-1| where p 𝑝 p is the optimal value returned by the solver (note that the matrix Y 𝑌 Y is sampled to have trace one).

$n$	time (s)	accuracy
5	2.37 s	1.143e-06
10	4.32 s	2.844e-06
15	9.56 s	4.732e-06
20	24.39 s	7.537e-06
25	77.02 s	9.195e-06
30	163.07 s	1.290e-05

Equations

$$\text{minimize}⟨c,x⟩\text{subject to}x\in L\cap {\mathbf{H}}_{+}^d$$

$$g(X)=U\mathrm{diag}(g({\lambda }_1),\dots ,g({\lambda }_n))U^{\ast }$$

$$g(\lambda X_1+(1-\lambda )X_2)⪰\lambda g(X_1)+(1-\lambda )g(X_2)$$

$$D(\rho \Vert \sigma ):=\mathrm{Tr}[\rho (\log \rho -\log \sigma )],$$

$$g(x)={\int }_t{f}_t(x)d\mu (t),$$

$$g(x)\approx \sum _{j=1}^m{w}_j{f}_{t_j}(x),$$

$$\log (x)={\int }_0^1\frac{x-1}{t(x-1)+1)}dt.$$

$$P_g(X,Y)=Y^{1/2}g\left(Y^{-1/2}X{Y}^{-1/2}\right)Y^{1/2}.$$

$$P_g\left(\lambda X_1+(1-\lambda )X_2,\lambda Y_1+(1-\lambda )Y_2\right)⪰\lambda P_g(X_1,Y_1)+(1-\lambda )P_g(X_2,Y_2)$$

$$\left\{(X,Y,T)\in {\mathbf{H}}_{++}^n\times {\mathbf{H}}_{++}^n\to {\mathbf{H}}^n:P_g(X,Y)⪰T\right\}.$$

$$D_{\text{op}}(X\Vert Y):=-X^{1/2}\log \left(X^{-1/2}Y{X}^{-1/2}\right)X^{1/2}.$$

$$\log (x)={\int }_1^x\frac{ds}{s}={\int }_0^1{f}_t(x)dt\text{ where }{f}_t(x)=\frac{x-1}{t(x-1)+1}.$$

$$f_t(x)\ge \tau ⟺\left[\begin{array}{cc}x-1-\tau & -\sqrt{t}\tau \\ -\sqrt{t}\tau & 1-t\tau \end{array}\right]⪰0.$$

$$\log (x)\approx \sum _{j=1}^m{w}_j{f}_{t_j}(x),$$

$$r_m(x):=\sum _{j=1}^m{w}_j{f}_{t_j}(x)=\sum _{j=1}^m{w}_j\frac{x-1}{t_j(x-1)+1}.$$

$$\log (x)=\frac{1}{h}\log (x^h).$$

$$r_{m,k}(x)=2^k{r}_m(x^{1/2^k}).$$

$$|r_{m,k}(x)-\log (x)|\le 2^k|\sqrt{\kappa }-\sqrt{{\kappa }^{-1}}{|}^2{\left(\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\right)}^{2m-1}\asymp 4\cdot 4^{-m(k+2)}\log (x{)}^{2m+1}(k\to \infty )$$

$$K_{\text{re}}:=\mathrm{cl}\left\{(x,y,\tau )\in {\mathbb{R}}_{++}\times {\mathbb{R}}_{++}\times \mathbb{R}:x\log (x/y)\le \tau \right\}.$$

$$K_{m,k}:=\{(x,y,t)\in {\mathbb{R}}_{++}^2\times \mathbb{R}:x{r}_{m,k}(x/y)\le t\}.$$

$$x{r}_{m,k}(x/y)\le x(\log (x/y)+ϵ)\le t+xϵ$$

$$f_t(X)⪰T\text{and}X\succ 0⟺\left[\begin{array}{cc}X-I& 0\\ 0& I\end{array}\right]-\left[\begin{array}{cc}T& \sqrt{t}T\\ \sqrt{t}T& tT\end{array}\right]⪰0\text{and}X\succ 0.$$

$$D_{\text{op}}(X\Vert Y):=-X^{1/2}\log \left(X^{-1/2}Y{X}^{-1/2}\right)X^{1/2}.$$

$$K_{\text{re}}^n=\mathrm{cl}\left\{(X,Y,T)\in {\mathbf{H}}_{++}^n\times {\mathbf{H}}_{++}^n\times {\mathbf{H}}^n:D_{\text{op}}(X\Vert Y)⪯T\right\}.$$

$$K_{m,k}^n=\left\{(X,Y,T)\in {\mathbf{H}}_{+}^n\times {\mathbf{H}}_{+}^n\times {\mathbf{H}}^n:-P_{r_{m,k}}(Y,X)⪯T\right\}$$

$$P_{r_{m,k}}(Y,X):=X^{1/2}{r}_{m,k}\left(X^{-1/2}Y{X}^{-1/2}\right)X^{1/2}.$$

$$\begin{array}[]{c}(X,Y,T)\in K_{m,k}^{n}\\ \rotatebox{90.0}{$\Longleftrightarrow$}\\ \exists T_{1},\dots,T_{m},Z_{0},\dots,Z_{k}\in\mathbf{H}^{n}\text{ s.t. }\left\{\begin{array}[]{ll}Z_{0}=Y,&\begin{bmatrix}Z_{i}&Z_{i+1}\\ Z_{i+1}&X\end{bmatrix}\succeq 0\;(i=0,\dots,k-1)\\[8.5359pt] \displaystyle\sum_{j=1}^{m}w_{j}T_{j}=-2^{-k}T,&\begin{bmatrix}Z_{k}-X-T_{j}&-\sqrt{t_{j}}T_{j}\\ -\sqrt{t_{j}}T_{j}&X-t_{j}T_{j}\end{bmatrix}\succeq 0\\ &\qquad\qquad\qquad(j=1,\dots,m)\end{array}\right.\end{array}$$

$$A{\#}_{h}B:=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^h{A}^{1/2}.$$

$$P_{r_{m,k}}(Y,X)=2^k{P}_{r_m}((X{\#}_{2^{-k}}Y),X).$$

$$Z_0=Y,Z_k=V\text{ and }\left[\begin{array}{cc}{Z}_i& Z_{i+1}\\ Z_{i+1}& X\end{array}\right]⪰0(i=0,\dots ,k-1).$$

Full text

1 Introduction

Semidefinite optimization problems are convex optimization problems that take the form

[TABLE]

where $\mathbf{H}^{d}_{+}$ is the cone of $d\times d$ Hermitian positive semidefinite matrices, and $L\subseteq\mathbf{H}^{d}$ is an affine subspace of $d\times d$ Hermitian matrices (thought of as a real vector space). A convex function $f$ is said to have a semidefinite representation of size $d$ if its epigraph $\{(x,t):f(x)\leq t\}$ can be expressed in the form $\pi(L\cap\mathbf{H}_{+}^{d})$ where $\pi$ is a linear map. The existence of such representations for many convex functions [BTN01a] explains the importance of semidefinite programming as a class of convex optimization problems. Understanding which convex sets and functions do and do not have small semidefinite descriptions has been a focus of considerable recent research effort in real algebraic geometry, optimization, and theoretical computer science (see, e.g., [BPT13]).

One fundamental limitation is that the feasible regions of semidefinite optimization problems are necessarily semialgebraic sets, i.e., they can be expressed as finite unions of sets defined by polynomial inequalities. As such, we cannot hope to exactly model non-semialgebraic convex sets and functions, such as the logarithm, using semidefinite programming. This leads us to consider the problem of understanding which general convex sets and functions can be approximated with high accuracy by sets with small semidefinite representations.

Semidefinite approximations

One starting point is to consider the approximation of univariate convex or concave functions from the point of view of semidefinite optimization. How well can we approximate a given univariate concave function with a function that is not just concave, but also has a semidefinite representation of a given size? This is distinct from questions in classical approximation theory, both due to its emphasis on preserving concavity, and also because the complexity of the approximating function is defined in terms of the size of a semidefinite description, rather than the degree of a polynomial or rational approximation. A key motivation for the study of univariate approximation theory is its relevance for computing matrix functions [Hig08, Tre13]. If $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ then the corresponding matrix function can be defined for positive definite Hermitian matrices $\mathbf{H}_{++}^{n}$ by

[TABLE]

where $X=U\operatorname{diag}(\lambda_{1},\ldots,\lambda_{n})U^{*}$ is an eigendecomposition of $X$. To generalize our semidefinite approximation point of view to matrix functions, we focus on functions that have a natural dimension-free concavity property known as operator concavity. A function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is operator concave if the corresponding matrix function satisfies Jensen’s inequality in the positive semidefinite (Löwner) order, i.e.,

[TABLE]

for all $n$, all $X_{1},X_{2}\in\mathbf{H}^{n}_{++}$ and all $\lambda\in[0,1]$. Associated with any operator concave function $g$ and a positive integer $n$ is a convex set $\{(X,T)\in\mathbf{H}_{++}^{n}\times\mathbf{H}^{n}:g(X)\succeq T\}$, the matrix hypograph of $g$. A good introduction to operator concave functions is [Car10].

Among the most familiar and important operator concave functions is the logarithm. The operator concavity of the logarithm has remarkable consequences. For example, it can be used to establish joint convexity of the (Umegaki) quantum relative entropy function,

[TABLE]

using an appropriate generalization of the perspective transform (the noncommutative perspective, to be defined later). The function $D$ plays a fundamental role in quantum information theory, and its joint convexity was first established by Lieb and Ruskai [LR73] building on an earlier result of Lieb [Lie73].

Contributions

In this paper we develop techniques to construct accurate approximations, with small semidefinite descriptions, for the matrix logarithm. A key motivation for doing so is that using this basic building block, we can approximate other important convex and concave functions arising in quantum information, such as the quantum relative entropy. The same basic principles we use to approximate the matrix logarithm apply in greater generality. Our methods partly generalize to yield high accuracy semidefinite approximations for functions that are operator monotone and operator concave, as well as their matrix analogues. Furthermore, the full power of our approximation methods for the matrix logarithm extend to operator concave functions that satisfy functional equations of a particular form. As examples in this direction we show how to obtain semidefinite approximations of the logarithmic mean, and the arithmetic-geometric mean of Gauss. We have implemented our constructions in the MATLAB-based modeling language CVX and they are available online on the website:

[TABLE]

Table 1 shows some of the functions implemented in the package.

Our functions can be combined with existing functions in CVX to solve problems involving a mixture of constraints modeled with the (operator) relative entropy cone, linear inequalities, and second-order and semidefinite cone constraints.

1.1 Key ideas

We now summarize the main ideas behind our approach to constructing semidefinite approximations and illustrate them with the central example of the paper, the logarithm.

Approximating integral representations via quadrature

The first main idea is to use integral representations of functions as the basis for approximation. In general, suppose a concave function $g$ has an integral representation of the form

[TABLE]

where $\mu$ is a positive measure and, for any fixed $t$, $f_{t}(x)$ is a semidefinite representable concave function of $x$. If we approximate the integral via a quadrature rule with positive weights (see Appendix A), we obtain an approximation of $g$ as

[TABLE]

which is again semidefinite representable. Integral representations of the form (3) are guaranteed to exist for certain operator concave functions by a result of Löwner. In the case of the logarithm, the integral representation is simply

[TABLE]

For fixed $t$, it turns out that the integrand is itself operator concave and its matrix hypograph has a semidefinite representation. Approximating the integral via a quadrature rule (such as Gaussian quadrature) we obtain an approximation of $\log$ that is operator concave and semidefinite representable.

Using functional equations to improve approximations

The logarithm also satisfies the functional equation $\log(x^{1/2})=\frac{1}{2}\log(x)$, allowing us to express $\log(x)$ in terms of the logarithm of $\sqrt{x}$. This is helpful because the square root brings points closer to $x=1$, where the approximations via quadrature are more accurate. Because the square root is also operator monotone, operator concave, and semidefinite representable, we can compose our rational approximations obtained via quadrature with this functional equation. Doing so we obtain improved approximations that still have all of these desirable properties.

This additional idea may seem specific to the logarithm. In fact, there are other operator monotone and operator concave functions obeying functional equations that relate the function at a point to the function value at a point closer to $x=1$. Moreover the functional equations have appropriate monotonicity and concavity properties, allowing us to use a similar strategy to obtain improved approximations. Functions defined as the limits of mean iterations, such as the arithmetic-geometric mean function of Gauss, have the appropriate properties to be approximated in this way.

Extending to bivariate matrix functions via perspectives

We can further extend our semidefinite approximations of matrix concave functions to certain bivariate matrix functions via a noncommutative notion of the perspective of a function. Given a function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$, its perspective transform is defined as $(x,y)\in\mathbb{R}^{2}_{++}\mapsto yg(x/y)$. It is a well-known result in convex analysis that if $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is concave then its perspective is also concave. The definition of the perspective transform extends to functions of positive definite matrices. Given a function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$, its noncommutative perspective is $P_{g}:\mathbf{H}_{++}^{n}\times\mathbf{H}_{++}^{n}\rightarrow\mathbf{H}^{n}$ defined by

[TABLE]

If $X$ and $Y$ are scalars, the noncommutative perspective coincides with the usual scalar definition of perspective transform. A remarkable property of the noncommutative perspective is that it is jointly concave in $(X,Y)$ whenever $g$ is operator concave, i.e.,

[TABLE]

for any $\lambda\in[0,1]$ and $X_{1},Y_{1},X_{2},Y_{2}\in\mathbf{H}^{n}_{++}$, see [Eff09, ENG11, EH14]. The semidefinite approximations we construct in this paper can be suitably homogenized to give semidefinite approximations of the noncommutative perspective, or more precisely of the associated hypograph cone:

[TABLE]

The noncommutative perspective of the negative logarithm function is known as operator relative entropy [FK89], which we denote by $D_{\text{op}}$:111We define $D_{\text{op}}$ as $D_{\text{op}}(X\|Y)=-P_{\log}(Y,X)$ to match the conventional order of arguments in information theory.

[TABLE]

The semidefinite approximations of the scalar logarithm function can be used to approximate $D_{\text{op}}$. In turn this allows us to get semidefinite approximations of the quantum relative entropy $D(\rho\|\sigma)=\operatorname{Tr}[\rho(\log\rho-\log\sigma)]$.

1.2 Related work

Computing the matrix logarithm

The problem of numerically computing the (matrix) logarithm has a long history in numerical analysis. Among the most successful methods is the so-called inverse scaling and squaring, or Briggs-Padé, method (see, e.g., [KL89, DMP96, AMH12]). This method uses the approximation $\log(X)\approx 2^{k}r_{m}(X^{1/2^{k}})$ where $r_{m}$ is the $m$th diagonal Padé approximant of $\log(x)$ at $x=1$, which turns out to be precisely the approximation we consider in this paper. The literature on computing the matrix logarithm via these methods does not seem to investigate the concavity properties of this approximation method. Our central observation is that this method for computing matrix logarithm “preserves” the concavity properties of logarithm, and can be modeled using semidefinite programming constraints. This in turn leads to efficient algorithms, via semidefinite programming, for problems much more complex than simply computing the matrix logarithm (see, e.g., [FF17]).

Other approximations

A simple approximation for logarithm is $\log(x)\approx\frac{1}{h}(x^{h}-1)$, where $0<h<1$, with equality when $h\rightarrow 0$. This can be seen as the combination of a Taylor linearization $\log(x)\approx x-1$ with the fact that $\log(x^{h})=\frac{1}{h}\log(x^{h})$. In previous work by the first two authors [FS16], it was shown that this approach can be used to get a semidefinite approximation of the matrix logarithm and the quantum relative entropy. In general, however, the quality of this approximation is relatively poor and is much less accurate than the rational approximations considered here. Another idea of approximating the scalar relative entropy cone via second-order cone programming is considered in unpublished work by Glineur [Gli09]. The approach taken by Glineur involves using an approximation for the logarithm via the arithmetic-geometric-mean iteration, and then giving an approximation of a convex cone related to the arithmetic-geometric-mean with convex quadratic inequalities.

Successive approximation

To make up for the poor approximation quality of $\log(x)\approx\frac{1}{h}(x^{h}-1)$, one method is to successively refine the linearization point and use, more generally, $\log(x)\approx\log(a)+\frac{1}{h}((x/a)^{h}-1)$. This is the approach taken by CVX [GB14]. It requires the solution of multiple second-order cone programs to update the linearization point. One drawback of this approach, however, is that it does not generalize to matrices, since there is no natural analogue of the identity $\log(ax)=\log(a)+\log(x)$ for matrices.

Approximating second-order cone programs with linear programs

The most prominent example of approximating a family of conic optimization problems with another, is the work of Ben-Tal and Nemirovski [BTN01b], giving a systematic method to approximate any second-order cone program with a linear program. The number of linear inequalities in the approximating linear programs of Ben-Tal and Nemirovski grow logarithmically with $1/\epsilon$ where $\epsilon$ is a notion of approximation quality. The fundamental construction underlying this approximation is a description of the regular $2^{n}$-gon in the plane as the projection of a higher-dimensional polyhedron with $2n$ facets. Using this technique Ben-Tal and Nemirovski give a polyhedral approximation to the exponential cone, via first constructing a second-order cone based approximation to the exponential cone [BTN01b, Example 4]. This approximation is based on a degree four truncation of the Taylor series for $\exp(2^{-k}x)$. Unlike our approximations, this approach works with the exponential, which is not operator convex and so does not generalize to matrices.

1.3 Outline

To make the presentation as accessible as possible, we focus first on the case of the logarithm function (Sections 2 and 3) before explaining the general approach for operator concave functions (Section 4). In Section 2 we describe the basic ideas behind our approximations, focusing on the scalar logarithm and the relative entropy cone. In Section 3 we state and prove our main result (Theorem 3), giving an explicit family of semidefinite approximations to the operator relative entropy. We conclude the section by giving semidefinite approximations of the epigraph of the quantum relative entropy function. In Section 4 we explain how our approach can be used to approximate other operator concave functions. In Section 5 we present some numerical experiments to test the accuracy of our approximations and give comparison with the successive approximation method of CVX. Finally we conclude in Section 6.

2 Approximating logarithm

In this section we describe the main ingredients for our semidefinite approximations of the logarithm. For simplicity we restrict ourselves, here, to the case of scalar logarithm. Nevertheles, our construction remains valid for matrices—we explain this is in more detail in the following section. Our approximation of the logarithm function relies on the following ingredients: an integral representation of $\log$, Gaussian quadrature, and the following functional relation satisfied by $\log$: $\log(x)=\frac{1}{h}\log(x^{h})$.

Integral representation

We start with the following integral representation of the logarithm function

[TABLE]

Here, the second equality comes from the change of variable $s=t(x-1)+1$. A key property of this integral representation is that for any fixed $t\in[0,1]$, the function $x\mapsto f_{t}(x)$ is concave. (The representation (6) thus establishes the concavity of $\log$ in a way that generalizes nicely to the setting of matrix functions.) One can easily show that the function $x\mapsto f_{t}(x)$ is semidefinite representable:

[TABLE]

Gaussian quadrature

To obtain an approximation of $\log$ that retains concavity, we discretize the integral (6) using Gaussian quadrature (see Appendix A for more information about Gaussian quadrature). This gives an approximation of the form

[TABLE]

where $t_{j}\in[0,1]$ are the quadrature nodes, and $w_{j}>0$ are the quadrature weights. We denote by $r_{m}(x)$ the right-hand side of (8), a rational function whose numerator and denominator have degree $m$:

[TABLE]

The key property of $r_{m}$ is that it is concave and semidefinite representable: this is because it is a nonnegative combination of functions that are each semidefinite representable (see (7)). It is also interesting to note that the function $r_{m}$ coincides precisely with the Padé approximant of $\log$ of type $(m,m)$: in particular $r_{m}$ agrees with the first $2m+1$ Taylor coefficients of the logarithm function. This has in fact been already observed, e.g., in [DMP96, Theorem 4.3] (see also Appendix B for a proof that works for a more general class of functions).

Exponentiation

The approximation (8) is best around $x=1$. A common technique to get good quality approximations when $x$ is farther away from $1$ is to exploit the following important property of the logarithm function:

[TABLE]

Note that when $0<h<1$, $x^{h}$ is closer to 1 than $x$ is, and thus the rational approximation (8) is of better quality at $x^{h}$ than at $x$. Taking $h$ of the form $h=1/2^{k}$ we define:

[TABLE]

The approximation $r_{m,k}$ should be understood as a composition of two steps for a given $x$: (1) take the $2^{k}$th root of $x$ to bring it closer to 1; and (2) apply the approximation $r_{m}$ and scale back by $2^{k}$ accordingly. One can show that $r_{m,k}$ is concave and semidefinite representable: indeed it is known that power functions of the form $x\mapsto x^{1/2^{k}}$ are concave and semidefinite representable (in fact second-order cone representable), see [BTN01a]. Since the function $r_{m}$ is concave, semidefinite representable, and monotone it easily follows that $r_{m,k}$ is concave and semidefinite representable. An explicit semidefinite representation appears as a special case of Theorem 3 in Section 3.

Error bounds

One can derive bounds on the error between $r_{m,k}$ and $\log$. Since $r_{m}$ is defined in terms of Gaussian quadrature applied to the rational function $f_{t}(x)$, such error bounds can be derived by studying the Chebyshev coefficients of $t\mapsto f_{t}(x)$. In fact these can be computed exactly and lead to the following error bounds.

Proposition 1**.**

Let $r_{m,k}$ be the function defined in (11). Then for any $x>0$ we have

[TABLE]

where $\kappa=x^{1/2^{k}}$.

Proof.

See Appendix B.2.2. ∎

By making appropriate choices of $m$ and $k$ in Proposition 1, we obtain a result showing how the size of our representation grows as the approximation quality improves.

Theorem 1**.**

For any (fixed) $a>1$ and any $\epsilon>0$, there exists a function $r$ such that $|r(x)-\log(x)|\leq\epsilon$ for all $x\in[1/a,a]$, and $r$ has a semidefinite representation of size $O(\sqrt{\log_{e}(1/\epsilon)})$.

Proof.

See Appendix B.2.2. ∎

The main point here is that it is the combination of Padé approximants with successive square rooting that allows us to get a rate of $O(\sqrt{\log(1/\epsilon)})$. Using either technique individually gives us a rate of $O(\log(1/\epsilon))$. Figure 1 shows the error $|r_{m,k}(x)-\log(x)|$ for different choices of $(m,k)$.

The (scalar) relative entropy cone

The relative entropy is defined as the perspective function of the negative logarithm: $(x,y)\in\mathbb{R}_{++}\times\mathbb{R}_{++}\mapsto x\log(x/y)$. The epigraph of this function is known as the relative entropy cone:

[TABLE]

Using the perspective of $r_{m,k}$ one can obtain a semidefinite approximation of $K_{\textup{re}}$. Let

[TABLE]

The following theorem gives an approximation error for the cone $K_{m,k}$.

Theorem 2** (Approximation error for $K_{\textup{re}}$).**

Let $a>1$ and $\epsilon>0$. Then there exist $m$ and $k$ with $m+k=O(\sqrt{\log_{e}(1/\epsilon)})$ such that:

if $0<a^{-1}y\leq x\leq ay$ and $(x,y,t)\in K_{\textup{re}}$ then $(x,y,t+x\epsilon)\in K_{m,k}$ 2. 2.

if $0<a^{-1}y\leq x\leq ay$ and $(x,y,t)\in K_{m,k}$ then $(x,y,t+x\epsilon)\in K_{\textup{re}}$.

Proof.

The proof is straightforward using Theorem 1. Theorem 1 says that there exist $(m,k)$ with $m+k=O(\sqrt{\log_{e}(1/\epsilon)})$ such that $|r_{m,k}(x)-\log(x)|<\epsilon$ on $[a^{-1},a]$. Now, if $(x,y,t)\in K_{\textup{re}}$ this means that $x\log(x/y)\leq t$. Since $x/y\in[a^{-1},a]$, we get that

[TABLE]

which means that $(x,y,t+x\epsilon)\in K_{m,k}$. The other direction is similar. ∎

A semidefinite representation of $K_{m,k}$ appears as the case $n=1$ of Theorem 3 to follow. Note that, in this scalar case, the approximation involves only $2\times 2$ linear matrix inequalities and thus can be formulated using second-order cone programming.

3 Operator concavity, noncommutative perspectives and the operator relative entropy cone

The main goal of this section is to show that the ideas presented in the previous section are still valid when working with matrices. The main result of this section (and of the paper) is Theorem 3, which gives an explicit semidefinite programming approximation of the operator relative entropy cone, a matrix generalization of the relative entropy cone.

We begin by showing that the approximation $r_{m,k}$, defined in (11), is operator concave, just like the logarithm function. We then show how to use the noncommutative perspective of $r_{m,k}$ to approximate the operator relative entropy. This leads to our explicit semidefinite approximation of the operator relative entropy cone. We then show how this can be used to approximate the quantum relative entropy.

3.1 Operator concavity of logarithm and its approximation

We have already mentioned in the introduction that the logarithm function is operator concave. The next proposition will allow us to show this, as well as the operator concavity of the rational function $r_{m}$ that we considered in the previous section.

Proposition 2**.**

For $t\in[0,1]$ let $f_{t}$ be the rational function defined in (6). Then $f_{t}$ is operator concave. In fact we have the following semidefinite representation of its matrix hypograph:

[TABLE]

Proof.

The fact that $f_{t}$ is operator concave will follow directly once we establish (12), since it will show that the matrix hypograph of $f_{t}$ is a convex set. The proof of (12) is based on Schur complements, and is given in Appendix C. ∎

We can now directly see that $\log$ is operator concave, since it is a nonnegative (integral) combination of the $f_{t}$. The same is also true for $r_{m}$, since it is defined as a finite nonnegative combination of the $f_{t}$. Since $f_{t}$ is semidefinite representable, we can also get a semidefinite representation of the matrix hypograph of $r_{m}$, i.e., $\{(Y,U):r_{m}(Y)\succeq U\}$. This is the special case of Theorem 3 with $k=0,U=-T$ and $X=I$.

Operator concavity of $r_{m,k}$

In Section 2 we saw that one can get an improved approximation of $\log$ by considering $r_{m,k}(x):=2^{k}r_{m}(x^{1/2^{k}})$. We now show that $r_{m,k}$ is also operator concave. The argument directly generalizes the proof that $r_{m,k}$ is concave (in the usual sense). For the generalization we need the notion of operator monotonicity. A function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is called operator monotone if whenever $X\succeq Y$ then $g(X)\succeq g(Y)$, where $X,Y\in\mathbf{H}^{n}_{++}$ for any $n$.

Proposition 3**.**

The function $r_{m,k}$ is operator concave.

Proof.

One can show that the functions $f_{t}$, for each fixed $t\in(0,1)$, are operator monotone in addition to being operator concave: this follows from the fact that $X\succeq Y\succ 0\implies X^{-1}\preceq Y^{-1}$. Since $r_{m}$ is a nonnegative combination of the $f_{t}$ it is also operator monotone and operator concave. It is well-known that the power functions $x\mapsto x^{1/2^{k}}$ are operator concave, see e.g., [Car10]. Finally it is not hard to show that the composition of an operator concave and monotone function, with an operator concave function, yields an operator concave function. Thus this proves operator concavity of $r_{m,k}$. ∎

We will see, by setting $X=I$ and $U=-T$ in Theorem 3 to follow, how to get an explicit semidefinite representation of the matrix hypograph, $\{(Y,U):r_{m,k}(Y)\succeq U\}$, of $r_{m,k}$.

3.2 Approximating the operator relative entropy cone

Recall that the operator relative entropy is the noncommutative perspective of the negative logarithm function:

[TABLE]

We know that $D_{\text{op}}$ is jointly matrix concave in $(X,Y)$. In particular, this means that the epigraph cone associated to $D_{\text{op}}$ is a convex cone:

[TABLE]

We saw, in Proposition 3, that $r_{m,k}$ is operator concave. It thus follows that the noncommutative perspective of $r_{m,k}$ is jointly concave. Our approximation of the cone $K_{\textup{re}}^{n}$ will be the epigraph cone of $-P_{r_{m,k}}$, the noncommutative perspective of $-r_{m,k}$. We will denote this cone by $K_{m,k}^{n}$:

[TABLE]

where

[TABLE]

The next theorem, which is the main result of this paper, gives an explicit semidefinite representation of the cone (15).

Theorem 3** (Main: semidefinite approximation of $K_{\textup{re}}^{n}$).**

The cone $K_{m,k}^{n}$ defined in (15) has the following semidefinite description:

[TABLE]

where $w_{j}$ and $t_{j}$ ($j=1,\dots,m$) are the weights and nodes for the $m$-point Gauss-Legendre quadrature on the interval $[0,1]$.

Proof.

To prove this theorem we need the notion of weighted matrix geometric mean. For $0<h<1$, the $h$-weighted matrix geometric mean of $A,B\succ 0$ is denoted $A\#_{h}B$ and defined by:

[TABLE]

Note that $A\#_{h}B$ is the noncommutative perspective of the power function $x\mapsto x^{h}$. The weighted matrix geometric mean is operator concave in $(A,B)$ and semidefinite representable. Semidefinite representations of $A\#_{h}B$ for any rational $h$ are shown in [Sag13, FS16].

Recall that $P_{r_{m}}$ is the noncommutative perspective of $r_{m}$. Since $r_{m,k}(X)=2^{k}r_{m}(X^{1/2^{k}})$, it is not difficult to verify that the noncommutative perspective $P_{r_{m,k}}$ of $r_{m,k}$ can be expressed as:

[TABLE]

The semidefinite representation (16) then follows from the following three facts:

Semidefinite representation of weighted matrix geometric means: For any $X,Y\succ 0$ and $V\in\mathbf{H}^{n}$ and $k\geq 1$ we have $X\#_{2^{-k}}Y\succeq V$ if and only if there exist $Z_{0},\dots,Z_{k}\in\mathbf{H}^{n}$ that satisfy:

[TABLE]

This is the case $h=1/2^{k}$ of the semidefinite representation that appears in [FS16]. This construction hinges on the fact that $X\#_{2^{-k}}Y$ can be expressed in terms of $k$ nested geometric means as $X\#_{1/2}(X\#_{1/2}(\dots(X\#_{1/2}Y)))$, the fact that

[TABLE]

and operator monotonicity of the geometric mean with respect to its arguments. 2. 2.

Semidefinite representation of $P_{r_{m}}$: For any $V,X\succ 0$ and $T\in\mathbf{H}^{n}$ we have $P_{r_{m}}(V,X)\succeq T$ if and only if there exist $T_{1},\dots,T_{m}$ that satisfy:

[TABLE]

This follows directly from the semidefinite representation of $P_{f_{t}}$ given in Proposition 8 (Appendix C) and the fact that $r_{m}=\sum_{j=1}^{m}w_{j}f_{t_{j}}$. 3. 3.

$P_{r_{m}}$ is monotone in its first argument. This easily follows from the monotonicity of $r_{m}$.

Combining these three ingredients, and using the expression of $P_{r_{m,k}}$ in Equation (18), yields the desired semidefinite representation (16). ∎

3.3 Quantum relative entropy

In this section, we see how to use the results from the previous section to approximate the (Umegaki) quantum relative entropy function, defined by

[TABLE]

The next proposition, which appears in [Tro15], shows how to express the epigraph of $D$ using the operator relative entropy cone (defined in Equation (14)).

Proposition 4** ([Tro15, Section 8.8]).**

Let $D$ be the relative entropy function (20) and $D_{\text{op}}$ be the operator relative entropy (13). Then for any $A,B\succ 0$ we have

[TABLE]

where $\phi$ is the unique linear map from $\mathbb{C}^{n^{2}\times n^{2}}$ to $\mathbb{C}$ that satisfies $\phi(X\otimes Y)=\operatorname{Tr}[XY^{T}]$, and $\bar{B}$ is the entrywise complex conjugate of $B$.

Proof.

We reproduce the proof in [Tro15, Section 8.8]. Observe that $A\otimes I$ and $I\otimes\bar{B}$ commute and, as such, $D_{\text{op}}(A\otimes I\|I\otimes\bar{B})=(A\otimes I)(\log(A\otimes I)-\log(I\otimes\bar{B}))$. Using the fact that $\log(X\otimes Y)=(\log X)\otimes I+I\otimes(\log Y)$, the previous equation simplifies to $D_{\text{op}}(A\otimes I\|I\otimes\bar{B})=(A\log A)\otimes I-A\otimes(\log\bar{B})$. Now, using the fact $\phi(X\otimes Y)=\operatorname{Tr}[XY^{T}]$, we immediately see that (21) holds. ∎

The previous proposition allows us to express the epigraph of the quantum relative entropy function (20) in terms of the operator relative entropy cone. This is the object of the next statement.

Corollary 1**.**

For any $A,B\succ 0$ and $\tau\in\mathbb{R}$ we have:

[TABLE]

Proof.

Straightforward from (21), and the fact that $X\preceq Y$ implies $\phi(X)\leq\phi(Y)$ (see Remark 1 below). ∎

One can then get a semidefinite approximation of the constraint $D(A\|B)\leq\tau$ by using the approximation given in (15) of the cone $K_{\textup{re}}^{n^{2}}$ and plugging it in (22). Note that the semidefinite approximation we thus get uses blocks of size $2n^{2}\times 2n^{2}$, because of the tensor product construction of Equation (21).

Remark 1**.**

Note that the linear map $\phi$ in Proposition 4 is given by $\phi(Z)=w^{*}Zw$ for $Z\in\mathbb{C}^{n^{2}\times n^{2}}$, where $w\in\mathbb{C}^{n^{2}}$ is the vector obtained by stacking the columns of the $n\times n$ identity matrix. It follows that $\phi$ is a positive linear map, in the sense that if $Z\succeq 0$ then $\phi(Z)\geq 0$.

A smaller representation

One can exploit the special structure of the linear map $\phi$ in (21), to reduce the size of the semidefinite approximation of $D(A\|B)$ from having $m+k$ blocks of size $2n^{2}\times 2n^{2}$, to having $m$ blocks of size $(n^{2}+1)\times(n^{2}+1)$ and $k$ blocks of size $2n^{2}\times 2n^{2}$. The main idea for this reduction is to observe that the rational function $f_{t}$, which is the main building block of our approximations, can be expressed as a Schur complement, namely we have $tP_{f_{t}}(X,Y)=Y-Y(Y+t(X-Y))^{-1}Y$. From this observation, one can get the following representation for the hypograph $v^{*}P_{f_{t}}(X,Y)v$, where $v\in\mathbb{C}^{n}$:

[TABLE]

This representation clearly has size $(n+1)\times(n+1)$. Combining this with the fact that $\phi$ has the form $\phi[X]=w^{*}Xw$ (see Remark 1), allows us to reduce the semidefinite approximation of $D(A\|B)$.

4 Approximating operator concave functions

The approximations to the relative entropy cone developed in Section 3 used the facts that

the logarithm is an integral of (semidefinite representable) rational functions, which can be approximated via quadrature; and 2. 2.

the logarithm obeys the functional equation $\log(\sqrt{x})=\frac{1}{2}\log(x)$.

In this section we show how to generalize these ideas, allowing us to give semidefinite approximations for convex cones of the form

[TABLE]

for a range of operator concave functions $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$, where $P_{g}(X,Y)$ is the noncommutative perspective of $g$ defined in (4). In Section 4.1, we discuss functions that admit similar integral representations to the logarithm, which can be approximated via quadrature. In Section 4.2, we present examples of functions with perspectives $P_{g}$ that obey functional equations of the form $P_{g}\circ\Phi=P_{g}$ where $\Phi$ is a map with certain monotonicity properties, and use these to obtain smaller semidefinite approximations.

4.1 Approximations via Löwner’s theorem

A general class of functions that admit integral representations are operator monotone functions, of which the logarithm is a special case. Recall that these are functions $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ that satisfy $g(X)\preceq g(Y)$ whenever $X\preceq Y$ for $X,Y\in\mathbf{H}^{n}_{++}$ and any $n\geq 1$. The following theorem, due to Löwner, shows that any operator monotone function admits an integral representation in terms of the rational functions $f_{t}$ that we saw earlier (see Appendix D).

Theorem 4** (Löwner).**

If $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is a non-constant operator monotone function then there is a unique probability measure $\nu$ supported on $[0,1]$ such that

[TABLE]

where $f_{t}$ is the rational function defined in (6).

The logarithm function corresponds to the case where the measure $\nu$, in (24), is the Lebesgue measure on $[0,1]$. One corollary of Löwner’s theorem is that any operator monotone function on $\mathbb{R}_{++}$ is necessarily operator concave, since the $f_{t}$ are operator concave, as we already saw. Given an operator monotone function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ (which is necessarily also operator concave) one can apply Gaussian quadrature on (24) (with respect to the measure $\nu$) to obtain a rational approximation of $g$. If we use $m$ quadrature nodes, we denote the corresponding rational function $r_{m}$. In Appendix B, we establish an error bound on the resulting approximation, which allows us to prove the following general theorem on semidefinite approximations of operator monotone functions.

Theorem 5** (Semidefinite approximation of operator monotone functions).**

Let $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ be an operator monotone (and hence operator concave) function and let $a>1$. Then for any $\epsilon>0$ there is a rational function $r$ such that $|r(x)-g(x)|\leq\epsilon$ for all $x\in[1/a,a]$, and $r$ has a semidefinite representation of size $O(\log(1/\epsilon))$.

Proof.

We show that $r=r_{m}$ has the desired properties. In Appendix B, Equation (42), we show that the error $|r_{m}(x)-g(x)|$ for $x\in[1/a,a]$ decays linearly in $m$, i.e., is $O(\rho^{m})$ for some constant $0<\rho<1$ depending on $a$. In other words if we take $m=O(\log(1/\epsilon))$ we get $|r_{m}(x)-g(x)|\leq\epsilon$ for all $x\in[1/a,a]$. Since each rational function $f_{t}$ has a semidefinite representation of size $2\times 2$ (see (7)) it follows that $r_{m}$ has a semidefinite representation of size $O(m)=O(\log(1/\epsilon))$ as a sum of $m$ such functions. ∎

The approximation $r$ we produce in Theorem 5 is also operator monotone and operator concave, and can be used to approximate the matrix hypograph of $g$, as well as its noncommutative perspective, just like for the logarithm function. The following result quantifies the error for the approximation of the cone $K_{g}^{n}$ (defined in (23)) we obtain this way.

Theorem 6** (Approximation error for $K_{g}^{n}$).**

Let $a>1$ and $\epsilon>0$ and let $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ be operator monotone (and hence operator concave). Then there exists $m$ with $m=O(\log_{e}(1/\epsilon))$ such that:

if $0\prec a^{-1}Y\preceq X\preceq aY$ and $(X,Y,T)\in K_{g}^{n}$ then $(X,Y,T+X\epsilon)\in K^{n}_{r_{m}}$ 2. 2.

if $0\prec a^{-1}Y\leq X\preceq aY$ and $(X,Y,T)\in K_{r_{m}}^{n}$ then $(X,Y,T+X\epsilon)\in K_{g}^{n}$.

Proof.

The proof is a straightforward matrix generalization of the proof of Theorem 2. We establish only the first statement, since the second is similar.

If $0\prec a^{-1}Y\preceq X\preceq aY$ then $a^{-1}I\preceq X^{-1/2}YX^{-1/2}\preceq aI$. By Theorem 5 there is $m=O(\log_{e}(1/\epsilon))$ such that $|r_{m}(x)-g(x)|\leq\epsilon$. Hence $-r_{m}(X^{-1/2}YX^{-1/2})+g(X^{-1/2}YX^{-1/2})\preceq\epsilon I$ and so, multiplying on the left and right by $X^{1/2}$, we see that $-P_{r_{m}}(X,Y)+P_{g}(X,Y)\preceq\epsilon X$. Since $(X,Y,T)\in K_{g}^{n}$, it follows that $-P_{g}(X,Y)\preceq T$ and so that $-P_{r_{m}}(X,Y)\preceq T+\epsilon X$. This shows that $(X,Y,T+\epsilon X)\in K_{r_{m}}^{n}$. ∎

Theorem 5 shows that, by just using Gaussian quadrature on (24), we can get semidefinite approximations of size $O(\log(1/\epsilon))$ for any operator monotone function $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$. In the next section we will see that if the function $g$ satisfies additional functional relations, then we can obtain approximations of order $O(\sqrt{\log(1/\epsilon)})$ or smaller. Before doing so, we consider certain positive-valued functions that will be useful later.

Positive-valued functions

In the special case when $g$ takes only positive values, one can prove (see Appendix D) an alternative integral representation, that has additional nice properties and takes the form

[TABLE]

Here, $\mu$ is a probability measure on $[0,1]$ and $f_{t}^{+}$ is the rational function $f_{t}^{+}(x)=((1-t)x^{-1}+t)^{-1}$. The main advantage of using this new integral representation instead of (24), is that the noncommutative perspective of $f_{t}^{+}$ is monotone with respect to both arguments, unlike $f_{t}$. This means that the perspective $P_{g}$ of any positive operator monotone function $g$, is monotone with respect to both arguments. Approximating $g$ by applying quadrature to (25) ensures this property is preserved.

Examples

The function $g(x)=x^{1/2}$ is known to be operator monotone and has the integral representation (25) with the measure $\mu$ given by the arcsine distribution:

[TABLE]

Another function known to be operator monotone is $g(x)=(x-1)/\log(x)$. In this case one can show that the measure $\mu$ in (25) is:

[TABLE]

More information about operator monotone functions and their integral representations can be found in the books by Bhatia [Bha09, Bha13].

4.2 Improved approximations via functional equations

The functional equation $\log(x^{1/2})=(1/2)\log(x)$ for the logarithm gives rise to a functional equation for the perspective, $P_{\log}(x,y)=y\log(x/y)$, of the logarithm. Indeed if we define

[TABLE]

In Section 2 we constructed rational approximations $r_{m}$ for the logarithm, and then improved the approximation quality by successive square-rooting, defining $r_{m,k}(x)=2^{k}r_{m}(x^{1/2^{k}})$. At the level of perspectives, we have that

[TABLE]

where $\Phi^{(k)}$ denotes the composition of $\Phi$ with itself $k$ times.

A similar approach is possible for operator monotone functions $g:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$, that satisfy a functional equation of the form $P_{g}\circ\Phi=P_{g}$, as long as $\Phi$ has certain monotonicity and contraction properties. In these cases, we can obtain semidefinite representable approximations to $g$ that have smaller descriptions, for a given approximation accuracy, than the approximations by rational functions given in Theorem 5. We make this precise in Theorem 7 to follow. For simplicity of notation, we work in the scalar setting, but our arguments all extend to the matrix setting.

Examples of operator monotone functions obeying a functional equation of the desired form come from the logarithmic mean and the arithmetic-geometric mean.

Logarithmic mean

In Section 4.1 we saw that the function $g(x)=\frac{x-1}{\log(x)}$ is operator monotone. Its perspective is the logarithmic mean:

[TABLE]

a function that arises naturally in problems of heat transfer, and in the Riemannian geometry of positive semidefinite matrices (see, e.g., [Bha09, Section 4.5]). If we define $\Phi(x,y)=((x+\sqrt{xy})/2,(y+\sqrt{xy})/2)$ then the logarithmic mean obeys the functional equation:

[TABLE]

The logarithmic mean also satisfies other functional equations that are closely related to Borchardt’s algorithm and variants [Car72] for computing the logarithm. These could also be used in the present context, but we focus on (28) for simplicity.

Arithmetic-geometric mean (AGM)

The arithmetic-geometric mean of a pair of positive scalars $x,y$, is defined as the common limit of the pair of (convergent) sequences $x_{0}=x$, $y_{0}=y$,

[TABLE]

This limit is denoted $\textup{AGM}(x,y)$, and is the perspective of the positive, operator monotone function, $g(x)=\textup{AGM}(x,1)$. Remarkably (see, e.g., [Cox04, Equation (1.7)]), the arithmetic-geometric mean is related to the complete elliptic integral of the first kind, $K(x)$, via

[TABLE]

Since it is defined as the limit of an iterative process, if $\Phi(x,y)=((x+y)/2,\sqrt{xy})$ then

[TABLE]

More examples can be obtained by considering operator monotone functions constructed via operator mean iterations, discussed, for instance, in [BP13].

4.2.1 Structure of approximations

Suppose $g:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$ is positive and operator monotone, and let $r_{m}^{+}$ be the rational, positive, operator monotone approximation to $g$ obtained by applying Gaussian quadrature (with respect to the measure $\mu$) to the integral representation (25). If, in addition, $P_{g}\circ\Phi=P_{g}$ for some map $\Phi:\mathbb{R}_{++}^{2}\rightarrow\mathbb{R}_{++}^{2}$, then we can define a two-parameter family of approximations by

[TABLE]

It makes sense to do this as long as $\Phi$ maps points ‘closer’ to the ray generated by $(1,1)$ (in a way made precise in Theorem 7, to follow), and the approximation $r_{m}^{+}$ of $g$ is accurate near $x=1$.

From now on we assume that $\Phi$ has the form

[TABLE]

where $h_{1},h_{2}:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$ are positive, operator monotone, functions. Observe that $\Phi$ has this form for the examples of the logarithmic mean and the arithmetic-geometric mean. If $\Phi$ has the form (30), then $P_{r_{m,k}}(x,y)$ (defined in (29)) is positive, jointly concave, and jointly monotone for all $k\geq 0$ and $m\geq 1$. In particular, these monotonicity and concavity properties ensure that the cones $K_{m,k}:=K_{r_{m,k}}$ can be (recursively) expressed as $K_{m,0}=K_{r_{m}^{+}}$ and

[TABLE]

for all $k\geq 1$. If the cones $K_{h_{1}}$ and $K_{h_{2}}$ associated with $h_{1}$ and $h_{2}$ have semidefinite descriptions of size $s_{1}$ and $s_{2}$ respectively, then $K_{m,k}$ has a semidefinite description of size $2m+k(s_{1}+s_{2})$.

4.2.2 Approximation error

The following result shows that if $\Phi$ has contraction and monotonicity properties, we can obtain smaller semidefinite approximations of nonnegative operator monotone functions $g$ satisfying a functional equation of the form $P_{g}\circ\Phi=P_{g}$. It allows us to get semidefinite approximations of size $O(\sqrt{\log(1/\epsilon)})$ if $\Phi$ contracts at a linear rate, and $O(\log\log(1/\epsilon))$ if $\Phi$ contracts quadratically, where $\epsilon$ is the approximation accuracy.

Theorem 7**.**

Let $g,h_{1},h_{2}:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$ be operator monotone (and hence operator concave) functions such that

[TABLE]

Suppose that $h_{1}$ and $h_{2}$ are semidefinite representable.

If there exists a constant $c>1$ such that

[TABLE]

then for any $a>1$ and any $\epsilon>0$ there is a function $r$ such that $|r(x)-g(x)|\leq\epsilon$ for all $x\in[1/a,a]$ and $r$ has a semidefinite representation of size $O(\sqrt{\log_{c}(1/\epsilon)})$.

If, in addition, there exists a constant $c_{0}>1$ such that

[TABLE]

then for any $a>1$ and any $\epsilon>0$ there is a function $r$ such that $|r(x)-g(x)|\leq\epsilon$ for all $x\in[1/a,a]$ and $r$ has a semidefinite representation of size $O(\log_{2}\log_{c_{0}}(1/\epsilon))$.

Proof.

We provide a proof in Appendix B.2.3. In each case we choose $r$ to be of the form $r_{m,k}(x)=P_{r_{m,k}}(x,1)$ (defined in (29)) for sufficiently large $m$ and $k$, and use the fact that $r_{m,k}$ has a semidefinite representation of size $O(m+k)$. ∎

Remark 2**.**

The condition (31) says that $d_{H}(\Phi(x,y),(1,1))\leq c^{-1}d_{H}((x,y),(1,1))$ where $d_{H}(\cdot,\cdot)$ is the Hilbert metric on rays of the cone $\mathbb{R}_{++}^{2}$ (see, e.g., [Bus73]). This is the precise sense in which $\Phi$ maps points ‘closer’ to the ray generated by $(1,1)$.

We now apply the theorem to the logarithmic mean and the arithmetic-geometric mean.

Logarithmic mean

In this case $g(x)=\frac{x-1}{\log(x)}$, $h_{1}(x)=(x+\sqrt{x})/2$, and $h_{2}(x)=(1+\sqrt{x})/2$. By a direct computation we see that

[TABLE]

Theorem 7 tells us that given $a>1$, there is a function $r:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$ with a semidefinite representation of size $O(\sqrt{\log_{2}(1/\epsilon)})$, such that $|r(x)-(x-1)/\log(x)|\leq\epsilon$ for all $x\in[1/a,a]$.

Arithmetic-Geometric mean

In this case $g(x)=\textup{AGM}(x,1)$, $h_{1}(x)=(x+1)/2$, and $h_{2}(x)=\sqrt{x}$. Then

[TABLE]

Furthermore, since $\log\cosh(z)\leq|z|$ for all $z$ and $\log\cosh(z)\leq z^{2}/2$ for all $z$, it follows that

[TABLE]

Theorem 7 tells us that given $a>1$, there is a function $r:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$ with a semidefinite representation of size $O(\log_{2}\log_{8}(1/\epsilon))$ such that $|r(x)-\textup{AGM}(x,1)|\leq\epsilon$ for all $x\in[1/a,a]$.

Remark 3**.**

As stated, both the construction of the functions $r_{m,k}$ in Section 4.2.1, and the statement of Theorem 7, are only valid when $g$ takes positive values. If $g$ is operator monotone but not positive-valued (as is the case for the logarithm), similar results apply if certain modifications are made. First, the rational functions $r_{m}$ (from Section 4.1) should be used in place of $r_{m}^{+}$ in (29). Second, we need the additional assumption that the second argument of $\Phi$ is linear (i.e., $h_{2}(x)$ is affine). This is required because $P_{r_{m}}$ is, in general, not monotone in its second argument.

5 Numerical experiments

We first evaluate our approximation method for the scalar relative entropy cone, and compare it with the successive approximation scheme of CVX, to solve maximum entropy problems and geometric programs. To assess the quality of the returned solutions, we use the solver Mosek [ApS15], which has a dedicated routine for entropy problems and geometric programming (mskenopt and mskgpopt respectively). Note, however, that this solver only deals with scalar problems, and has no facility for matrix problems involving quantum relative entropy, for instance. To evaluate our method for matrices, we test it on a variational formula for trace. More numerical experiments using CvxQuad related to problems in quantum information theory appear in [FF17].

5.1 Entropy problems

We consider optimization problems of the form

[TABLE]

and we compare the performance of our method with the successive approximation scheme implemented in CVX. Table 2 shows the results of the comparison for randomly generated data $A\in\mathbb{R}^{\ell\times n}$ and $b\in\mathbb{R}^{\ell}$ of different sizes. We use the solution returned by the built-in maximum entropy solver in Mosek (mskenopt) as “true solution” and we measure the quality of either approximation method (successive approximation or ours) via the gap between optimal values. We use the notation $p_{sa}$ and $p_{Pade}$ respectively for the optimal values returned by the successive approximation scheme and our method.

5.2 Geometric programming

We consider now geometric programming [BKVH07] problems of the form:

[TABLE]

where $x\in\mathbb{R}^{n}$ is the decision variable. For $a\in\mathbb{R}^{n}_{++}$ the notation $x^{a}$ indicates $x^{a}:=\prod_{i=1}^{n}x_{i}^{a_{i}}$. The coefficients $c_{j,k}$ are assumed to be positive. Such problems can be converted into conic problems over the relative entropy (exponential) cone using the change of variables $y_{i}=\log x_{i}$. The current version of CVX (CVX 2.1) uses the successive approximation technique to deal with such problems. Our method based on Padé approximations can also be used in this case to obtain accurate approximations. We note that the solver Mosek has a dedicated routine for geometric programming (mskgpopt).

Table 3 shows a comparison of our method with the successive approximation method for randomly generated instances of (34). The instances were generated using the mkgp script contained in the ggplab package available at https://stanford.edu/~boyd/ggplab/.

5.3 Variational formula for trace

We now evaluate our method for matrix functions. We consider the following variational expression for the trace function which appears in [Tro12, Lemma 6]. For any $Y\succ 0$

[TABLE]

where $D$ is the quantum relative entropy function (2). We generate random positive definite matrices $Y$ and compare the solution of the right-hand side of (35) with $\operatorname{Tr}[Y]$. The right-hand side of (35) can be implemented using the CVX code shown in Table 4. The results of running this piece of code using solver SDPT3 are shown in Table 4.

6 Discussion

Lower bounds

It would be interesting to know what is the smallest possible second-order cone program that can approximate logarithm to within a fixed $\epsilon>0$. To formalize this question, let $\mathcal{F}_{s}$ be the class of concave functions on $\mathbb{R}_{++}$ that admit a second-order cone representation of size at most $s$.

[TABLE]

Recall, from Theorem 1, that our construction yields $s(\epsilon)=O(\sqrt{\log(1/\epsilon)})$. This rate results from the combination of Padé approximation with successive square rooting. It would be interesting to produce lower bounds on $s(\epsilon)$.

More generally one can define a notion of $\epsilon$-approximate extension complexity of a concave function $g:[a,b]\rightarrow\mathbb{R}$ in a similar way as (36). Well-known results in classical approximation theory relate the approximation quality using polynomials and rational functions of given degree to the smoothness of $g$. A natural question is to understand what corresponding properties of a concave function make it more or less difficult to approximate using second-order programs. We have phrased the question here in terms of second-order cone representations for concreteness but the same question for linear programming and semidefinite programming can also be considered.

Smaller semidefinite approximations for quantum relative entropy

The approximations for the epigraph of the quantum relative entropy $D(A\|B)$ we constructed in Section 3.3 involve linear matrix inequalities of size $O(n^{2})$ (where $n$ is the size of the matrices $A,B$). Is it possible to obtain approximations, of similar quality, to the quantum relative entropy using linear matrix inequalities of size $O(n)$?

Self-concordant barriers for the operator relative entropy cone

A natural approach to conic optimization over the scalar relative entropy cone (or, equivalently, the exponential cone) is to use an interior point method that works directly with an efficiently computable self-concordant barrier for the cone (such as the barrier introduced by Nesterov [Nes06]). Examples of such solvers include the extension of ECOS [DCB13] to the exponential cone [Ser15], and the solver developed by Skajaa and Ye [SY15]. We are not aware, however, of any barrier for the operator relative entropy cone that is known to be efficiently computable and self-concordant. If we had such a barrier, it could be used directly to solve conic optimization problems over the operator relative entropy cone using interior point methods, as an alternative to the semidefinite approximation-based approaches developed in this paper.

Approximating other families of convex functions via quadrature

One of the basic ideas of this paper is that if we can express a convex (or concave) function as $g(x)=\int_{\alpha}^{\beta}K(x,t)\;d\mu(t)$, where $x\mapsto K(x,t)$ has a simple semidefinite representation for fixed $t$, then we can obtain a semidefinite approximation of $g$ by quadrature. Operator monotone functions on $\mathbb{R}_{++}$, such as the logarithm, are just one class of functions with such a representation. Other such families of functions include Stieltjes functions, and certain hypergeometric functions. For instance Stieltjes functions on $\mathbb{R}_{++}$ have the form $g(x)=\int_{0}^{\infty}\frac{1}{x+t}d\mu(t)$. Hypergeometric functions ${}_{2}F_{1}(a,b;c;x)$ for $x<1$ and $b,c>0$ have the form ${}_{2}F_{1}(a,b;c;x)=\frac{1}{B(b,c-b)}\int_{0}^{1}t^{b-1}(1-t)^{c-b-1}(1-xt)^{-a}\;dt$ where $B(\cdot,\cdot)$ is the beta function. We expect such integral representations to be helpful in the study of approximate extension complexity of functions.

Free semidefinite representation

The semidefinite representation given in this paper of the hypograph of $f_{t}$ (see (12)) is a “free linear matrix inequality” representation in the sense of [HKM17]. This is one reason why our representations also work for the noncommutative perspective of $f_{t}$. In fact one can show that if an operator concave function $f$ admits a free linear matrix inequality representation, then the noncommutative perspective of $f$ also has a free linear matrix inequality representation. An interesting question would be to understand the class of operator concave functions that admit a free LMI representation.

Appendix A Background on approximation theory

Gaussian quadrature

A quadrature rule is a method of approximating an integral with a weighted sum of evaluations of the integrand. A quadrature rule is determined by the evaluation points, called nodes, and the weights of the weighted sum. Given a measure $\nu$ supported on $[-1,1]$, a quadrature rule gives an approximation of the form

[TABLE]

where the $t_{j}\in[-1,1]$ are the nodes and the $w_{j}$ are the weights. A Gaussian quadrature is a choice of nodes $t_{1},\ldots,t_{m}\in(-1,1)$ and positive weights $w_{1},\ldots,w_{m}$ that integrates all polynomials of degree at most $2m-1$ exactly. For example, when $\nu$ is the uniform measure on $[-1,1]$, such a quadrature rule is known as Gauss-Legendre quadrature, and the nodes and weights can be computed for example by an eigenvalue decomposition of the associated Jacobi matrix, see e.g., [Tre08, Section 2].

Padé approximants

Padé approximants are approximations of a given univariate function, analytic at a point $x_{0}$, by rational functions. More precisely, the $(m,n)$-Padé approximant of $h$ at $x_{0}$ is the rational function $p(x)/q(x)$ such that $p$ is a polynomial of degree $m$, $q$ is a polynomial of degree $n$, and the Taylor series expansion of the error at $x_{0}$ is of the form

[TABLE]

for real numbers $a_{k}$ and the largest possible positive integer222This last requirement is to ensure uniqueness. $s$. Expressed differently, $p$ and $q$ are chosen so that the Taylor series of $p(x)/q(x)$ at $x_{0}$ matches as many Taylor series coefficients of $h$ at $x_{0}$ as possible (and at least the first $m+n+1$ coefficients).

Appendix B Properties and error bounds of Gaussian quadrature-based approximations

Assume $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is a function with an integral representation

[TABLE]

where $\nu$ is a probability measure on $[0,1]$ and $f_{t}(x)=\frac{x-1}{1+t(x-1)}$. The case $g=\log$ corresponds to $\nu$ being the Lebesgue measure on $[0,1]$. In this appendix we show that the rational approximation obtained by applying Gaussian quadrature on (38) coincides with the Padé approximant of $g$ at $x=1$. We also derive error bounds on the quality of this rational approximation. Note that functions of the form (38) are precisely operator monotone functions, by Löwner’s theorem (see Section 4.1).

Let $r_{m}$ be the rational approximant obtained by using Gaussian quadrature on (38):

[TABLE]

where $w_{i}>0,t_{i}\in[0,1]$ are the Gaussian quadrature weights and nodes for the measure $\nu$.

B.1 Connection with Padé approximant

We first show that the function $r_{m}$ coincides with the Padé approximant of $g$ at $x=1$. The special case $g=\log$ was established in [DMP96, Theorem 4.3].

Proposition 5**.**

Assume $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ has the form (38) and let $r_{m}$ be the rational approximation obtained via Gaussian quadrature as in (39). Then $r_{m}$ is the $(m,m)$ Padé approximant of $g$ at $x=1$.

Proof.

First we note that $f_{t}(x)$ admits the following series expansion, valid for $|x-1|<\frac{1}{|t|}$:

[TABLE]

Let $\nu_{m}=\sum_{i=1}^{m}w_{i}\delta_{t_{i}}$ be the atomic measure on $[0,1]$ corresponding to Gaussian quadrature applied to $\nu$. By definition of Gaussian quadrature, $\nu_{m}$ matches all moments of $\nu$ up to degree $2m-1$, i.e., $\int_{0}^{1}p(t)\;d\nu(t)=\int_{0}^{1}p(t)\;d\nu_{m}(t)$ for all polynomials $p$ of degree at most $2m-1$. It thus follows that

[TABLE]

establishing that $r_{m}$ matches the first $2m$ Taylor series coefficients of $g$ at $x=1$. Since $r_{m}$ has numerator and denominator degree $m$, it is the $(m,m)$-Padé approximant of $g$ at $x=1$. ∎

B.2 Error bounds

In this section we derive an error bound on the approximation quality $|g(x)-r_{m}(x)|$. To do this we use standard methods as described, e.g., in [Tre13]. This error is essentially controlled by the decay of the Chebyshev coefficients of the integrand. For the rational functions $f_{t}$ one can compute these coefficients exactly.

B.2.1 Quadrature error bounds for operator monotone functions

To appeal to standard arguments, it is easiest to rewrite the integrals of interest over the interval $[-1,1]$ by the transformation $t\mapsto 1-2t$ mapping $[0,1]$ to $[-1,1]$. To this end, let

[TABLE]

Let $T_{k}(t)$ denote the $k$th Chebyshev polynomial. We start by explicitly computing the Chebyshev expansion of $\tilde{f}_{t}(x)$ for fixed $x$, i.e., we find the coefficients $a_{k}(x)$ of $\tilde{f}_{t}(x)=\sum_{k=0}^{\infty}a_{k}(x)T_{k}(t)$. To do this, we first define $h_{\rho}(t)=\frac{2}{(\rho+\rho^{-1})/2-t}$ and observe that with the substitution $\rho=\frac{\sqrt{x}-1}{\sqrt{x}+1}$ we have that $\tilde{f}_{t}(x)=h_{\rho}(t)$ and that $x>0$ if and only if $-1<\rho<1$. We can compute the Chebyshev expansion of $h_{\rho}(t)$ by observing that the generating function of Chebyshev polynomials is (see e.g., [Tre13, Exercise 3.14])

[TABLE]

It then follows that the Chebyshev expansion of $h_{\rho}(t)$ is

[TABLE]

Since $\frac{8}{\rho^{-1}-\rho}=2(\sqrt{x}-1/\sqrt{x})$, the Chebyshev expansion of $\tilde{f}_{t}(x)$ is

[TABLE]

We are now ready to state an error bound on the approximation quality $|g(x)-r_{m}(x)|$. Our arguments are standard, and follow closely the ideas described in [Tre13].

Proposition 6**.**

Let $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ be a function with an integral representation (38) and let $r_{m}$ be the rational approximation obtained by applying Gaussian quadrature as in (39). If $m\geq 1$ and $x>0$ then

[TABLE]

If $\nu$ is invariant under the map $t\mapsto 1-t$ (i.e., $g(x^{-1})=-g(x)$) then this can be improved to

[TABLE]

Finally, $r_{m}(x)\geq g(x)$ for all $0<x\leq 1$ and $r_{m}(x)\leq g(x)$ for all $x\geq 1$.

Proof.

Let $\tilde{\nu}$ be the measure on $[-1,1]$ obtained from $\nu$ by changing variables $t\in[0,1]\mapsto 1-2t\in[-1,1]$ so that $g(x)=g(0)+g^{\prime}(1)\int_{-1}^{1}\tilde{f_{t}}(x)d\tilde{\nu}(t)$. Let $\tilde{\nu}_{m}$ be the atomic measure supported on $m$ points obtained by applying Gaussian quadrature on $\nu$. Finally let the Chebyshev expansion of $\tilde{f}_{t}(x)$ be $\sum_{k=0}^{\infty}a_{k}(x)T_{k}(t)$. Since $\int_{-1}^{1}T_{k}(t)\;d\tilde{\nu}(t)=\int_{-1}^{1}T_{k}(t)\;d\tilde{\nu}_{m}(t)$ for $k\leq 2m-1$,

[TABLE]

For $k\geq 2$, we have that $a_{k}(x)=2(\sqrt{x}-1/\sqrt{x})\left(\frac{\sqrt{x}-1}{\sqrt{x}+1}\right)^{k}$ (see (41)). So using the fact that $\tilde{\nu}$ and $\tilde{\nu}_{m}$ are probability measures (when $m\geq 1$), together with the fact that $|T_{k}(t)|\leq 1$ for $t\in[-1,1]$, the triangle inequality gives

[TABLE]

If the measure $\nu$ is invariant under the map $t\mapsto-t$ then the same is true of $\nu_{m}$ (see, e.g., [MS14]). Since $\tilde{f}_{t}(x^{-1})=-\tilde{f}_{-t}(x)$ it follows that $r_{m}(x^{-1})=-r_{m}(x)$. Furthermore,

[TABLE]

because Chebyshev polynomials of odd degree are odd functions. In this case only the even Chebyshev coefficients contribute to the error bound so

[TABLE]

To establish inequalities between $r_{m}(x)$ and $g(x)$, we use an alternative formula for the error obtained by approximating an integral via Gaussian quadrature. Since $t\mapsto\tilde{f}_{t}(x)$ has derivatives of all orders, one can show (see, e.g., [SB02, Theorem 3.6.24]) that there exists $\tau\in[-1,1]$ and $\kappa\geq 0$ such that

[TABLE]

If $x\in(0,1)$ then $\frac{1+x}{1-x}-\tau<0$ for all $\tau\in[-1,1]$ and so $g(x)-r_{m}(x)<0$. If $x\in(1,\infty)$ then $\frac{1+x}{1-x}-\tau>0$ for all $\tau\in[-1,1]$ and so $g(x)-r_{m}(x)>0$. If $x=1$ then $g(x)=r_{m}(x)$. ∎

Very similar bounds hold for the error between $g:\mathbb{R}_{++}\rightarrow\mathbb{R}_{++}$, a positive operator monotone function, and $r_{m}^{+}$, the rational approximation obtained by applying Gaussian quadrature to the integral representation in (25). Indeed if $m\geq 1$ and $x>0$,

[TABLE]

We omit the proof, since it follows the same basic argument as the proof of (42), together with the observation that $f_{t}^{+}(x)=\frac{x}{x-1}f_{t}(x)$.

B.2.2 The special case of log: proofs of Proposition 1 and Theorem 1

Proof of Proposition 1.

The function $g(x)=\log(x)$ has an integral representation (38) where the measure $\nu$ is the uniform measure on $[0,1]$, which is invariant under the map $t\mapsto 1-t$. Proposition 6 tells us that, for any $x>0$,

[TABLE]

The error between $\log(x)=2^{k}\log(x^{1/2^{k}})$ and $r_{m,k}(x)=2^{k}r_{m}(x^{1/2^{k}})$ can be obtained by evaluating at $x^{1/2^{k}}$ and scaling by $2^{k}$ to obtain

[TABLE]

where $\kappa=x^{1/2^{k}}$. By using the fact that

[TABLE]

we can write this as a bound on relative error as

[TABLE]

Asymptotic behavior of (46): Since $\kappa=x^{1/2^{k}}=e^{2^{-k}\log(x)}$, we can rewrite the right-hand side of (46) as

[TABLE]

Since $\sinh^{2}(2x)\tanh^{2m-1}(x)=4x^{2m+1}+O(x^{2m+3})$, we have that

[TABLE]

∎

Proof of Theorem 1.

The function $r$ is chosen to be of the form $r_{m,k}$ for certain $m$ and $k$. In particular we can choose $k=k_{1}+k_{2}$, with $k_{1}=\lceil\log_{2}\log_{e}(a)\rceil+1$, $k_{2}$ being the smallest even integer larger than $\sqrt{\log_{2}(32\log_{e}(a)/\epsilon)}$, and with $m=k_{2}/2$. The function $r_{m,k}$ has a semidefinite representation of size $m+k$ (as a special case of Theorem 3 in Section 3), which is $O(\sqrt{\log_{e}(1/\epsilon)})$ for fixed $a$. It remains to establish the error bound. To do so, we first note that $x^{1/2^{k_{1}}}<1$ for all $x\in[1/a,a]$. Then, for all $x\in[1/a,a]$,

[TABLE]

Here, the second last equality holds because $\sinh(1/2)^{2}\leq 1$, $\tanh(x)\leq x$ for all $x\geq 0$, and $2^{k_{1}-1}\leq\log_{e}(a)$ (by our choice of $k_{1}$). The last inequality holds by our choice of $m$ and $k_{2}$. ∎

B.2.3 Proof of Theorem 7

Proof of Theorem 7.

The function $r$ is of the form $r_{m,k}$ defined by (29) for particular values of the parameters $m$ and $k$. Throughout the proof, for convenience of notation, let $(x_{k},y_{k})=\Phi^{(k)}(x,y)$ for $k\geq 0$. The error bound (44) of Appendix B.2 shows that for any $x,y>0$:

[TABLE]

We will show that if $\Phi$ has the linear contraction property (31) then the bound (47) decays like $O(c^{-k^{2}})$ for the choice of $m\approx k$ (that we make precise later). To establish this, we need to bound two terms: first, if (31) holds then $\log(x_{k}/y_{k})=O(c^{-k/2})$ and so the numerator in (47) converges like $O(c^{-km})$ as we want. The second term that we need to control is $\sqrt{x_{k}y_{k}}$ and one can show that this term grows at most linearly. This is proved in the following lemma:

Lemma 1**.**

There is a constant $b>0$ such that for any $x,y>0$ satisfying $a^{-1}\leq x/y\leq a$ we have

[TABLE]

where $(x_{k},y_{k})=\Phi^{(k)}(x,y)$.

Proof.

Since $h_{1}$ and $h_{2}$ are concave, they are each bounded above by their linear approximation at $x=1$. As such, $P_{h_{i}}(x,y)\leq h_{i}^{\prime}(1)(x-y)+h_{i}(1)y$ for all $x,y\in\mathbb{R}_{++}$ and $i=1,2$. Summing these two inequalities we see that

[TABLE]

Because $h_{1}$ and $h_{2}$ take positive values, $h_{1}(1)\geq h_{1}(1)-h_{1}^{\prime}(1)\geq 0$ and $h_{2}(1)\geq h_{2}(1)-h_{2}^{\prime}(1)\geq 0$. As such, if $b=\max\{h_{1}^{\prime}(1)+h_{2}^{\prime}(1),h_{1}(1)+h_{2}(1)-(h_{1}^{\prime}(1)+h_{2}^{\prime}(1))\}$, then $P_{h_{1}}(x,y)+P_{h_{2}}(x,y)\leq b(x+y)$ for all $x,y\in\mathbb{R}_{++}$. It then follows that $x_{k}+y_{k}\leq b^{k}(x+y)$ for all $x,y\in\mathbb{R}_{++}$ and so that

[TABLE]

as desired. ∎

Plugging (48) in (47) gives us, for any $a^{-1}\leq x/y\leq a$:

[TABLE]

Choose $k$ to be the smallest even integer satisfying $k\geq\max\left\{2\log_{c}\log(a),\sqrt{\log_{c}\left(\frac{8(g(1)-g(0))(1+a)}{3\epsilon}\right)}\right\}$ and $m$ to be the smallest integer satisfying $m\geq k\max\{1,\frac{\log(b)}{\log(16)}\}$. Note that both $m$ and $k$ are $O(\sqrt{\log_{c}(1/\epsilon)})$ when we treat $a$ and $b$ as constants. With these choices, and the assumption (31), we have that $b^{k}16^{-m}\leq 1$ and

[TABLE]

Using the inequality $|\tanh(z)|\leq|z|$ for all $z$, and setting $y=1$ in the error bound (49), we have that

[TABLE]

The size of the semidefinite representation of $r=r_{m,k}$ is $O(m+k)$, if we view the size of the semidefinite representations of $h_{1}$ and $h_{2}$ as being constant. Since $m,k\in O(\sqrt{\log_{c}(1/\epsilon)})$ it follows that the size of the semidefinite representation of $r$ is also $O(\sqrt{\log_{c}(1/\epsilon)})$.

In the case where the assumption (32) also holds, we choose $m$ (respectively $k$) to be the smallest integer (respectively even integer) satisfying

[TABLE]

Note that both $m$ and $k$ are $O(\log_{2}\log_{c_{0}}(1/\epsilon))$. With these choices, and the assumptions (31) and (32), we have that $c_{0}^{m}b^{k}16^{-m}\leq 1$ and $|\log(x_{k/2}/y_{k/2})|\leq c^{-k/2}\log(a)\leq 1$ and

[TABLE]

Using the inequality $|\tanh(z)|\leq|z|$ for all $z$, and putting $y=1$ in the error bound (49), we obtain

[TABLE]

∎

Appendix C Semidefinite description of $f_{t}$

In this section we establish the linear matrix inequality characterization of $f_{t}$ given in Proposition 2. We use the fact that if $t\in(0,1]$ then

[TABLE]

The characterization will follow from the following easy observation.

Proposition 7**.**

If $A+B\succ 0$ then

[TABLE]

Proof.

The proof follows by expressing the left-hand side of (51) using Schur complements, followed by a congruence transformation:

[TABLE]

∎

Proof of Proposition 7.

We need to show that

[TABLE]

The case $t=0$ can be easily verified to hold. We thus assume $0<t\leq 1$. Given the expression of $f_{t}$ in Equation (50) we simply apply (51) with $B=(I/t)$ and $A=X-I$. This shows that

[TABLE]

Applying a congruence transformation with the diagonal matrix $\operatorname{diag}(I,\sqrt{t}I)$ yields the desired linear matrix representation (52). ∎

We can also directly get, from Proposition 7, a semidefinite representation of the noncommutative perspective of $f_{t}$ defined by $P_{f_{t}}(X,Y)=Y^{1/2}f_{t}\left(Y^{-1/2}XY^{-1/2}\right)Y^{1/2}$.

Proposition 8**.**

If $t\in[0,1]$ then the perspective $P_{f_{t}}$ of $f_{t}$ is jointly matrix concave since

[TABLE]

Proof.

From the definition of $P_{f_{t}}$ and the expression (50) for $f_{t}$ it is easy to see that we have:

[TABLE]

The semidefinite representation (53) then follows easily by applying (51) with $B=Y/t$ and $A=X-Y$, followed by applying a congruence transformation with the diagonal matrix $\operatorname{diag}(1,\sqrt{t})$. ∎

Appendix D Integral representations of operator monotone functions

In this section we show how to obtain the integral representations (24) and (25) as a fairly easy reworking of the following result.

Theorem 8** ([HP82, Theorem 4.4]).**

If $h:(-1,1)\rightarrow\mathbb{R}$ is non-constant and operator monotone then there is a unique probability measure $\tilde{\nu}$ supported on $[-1,1]$ such that

[TABLE]

Suppose $g:\mathbb{R}_{++}\rightarrow\mathbb{R}$ is operator monotone. Then it is straightforward to check that $h:(-1,1)\rightarrow\mathbb{R}$ defined by $h(z)=g\left(\frac{1+z}{1-z}\right)$ is operator monotone and that $g(x)=h\left(\frac{x-1}{x+1}\right)$. By applying Theorem 8 to $h(z)$ and then evaluating at $z=\frac{x-1}{x+1}$, we obtain the integral representation

[TABLE]

Using the fact that $h(0)=g(1)$ and $h^{\prime}(0)=2g^{\prime}(1)$, and applying a linear change of coordinates to rewrite the integral over $[0,1]$, we see that there is a probability measure $\nu$ on $[0,1]$ such that

[TABLE]

This establishes (24). If, in addition, $g$ takes positive values, then $g(0):=\lim_{x\rightarrow 0}g(x)\geq 0$. Hence

[TABLE]

so we can define a probability measure supported on $[0,1]$ by $d\mu(t)=\frac{g^{\prime}(1)}{g(1)-g(0)}\left(\frac{1}{1-t}\right)\;d\nu(t)$. Then using the fact that $f_{t}(x)=\frac{1}{1-t}\left[f^{+}_{t}(x)-1\right]$ we immediately obtain, from (55), the representation

[TABLE]

This establishes (25).