Matrix Functions and Scaling Squaring

For an integer $m\ge 0$, put $a_{m,j}:=\frac{(2m-j)!m!}{(2m)!j!(m-j)!}(0\le j\le m),p_m(z):={\sum }_{j=0}^m{a}_{m,j}{z}^j,q_m(z):=p_m(-z)\mathrm{.}$ The diagonal Padé rational function of order $m$ for the scalar exponential The complex exponential by its power series is $r_m(z):=\frac{p_m(z)}{q_m(z)}$ at every $z\in \mathbb{C}$ for which $q_m(z)\mathrm{\ne }0$. The pair is normalized by $p_m(0)=q_m(0)=1$. This definition does not assign a value at a zero of the denominator.

Let $s\ge 0$, let ${\lambda }_1,\dots ,{\lambda }_s\in \mathbb{C}$ be distinct, let $m_1,\dots ,m_s$ be positive integers, and put $N={\sum }_{i=1}^s{m}_i$. Let $V_N=\mathrm{span}\{1,z,\dots ,z^{N-1}\}$, with the empty span $V_0=\{0\}$. For arbitrary data $c_{i,k}\in \mathbb{C}$ with $1\le i\le s$ and $0\le k<m_i$, there is exactly one polynomial $p\in V_N$ satisfying $p^{(k)}({\lambda }_i)=c_{i,k}(1\le i\le s,0\le k<m_i)\mathrm{.}$

Let $\mathrm{\parallel }\cdot \mathrm{\parallel }$ be a submultiplicative norm on square matrices. Suppose $X_{j+1}=X_j^2$ and ${\widehat{X}}_{j+1}={\widehat{X}}_j^2+E_j$ for matrices of the same size. If $e_j=\mathrm{\parallel }{\widehat{X}}_j-X_j\mathrm{\parallel }$, then $e_{j+1}\le (2\mathrm{\parallel }{X}_j\mathrm{\parallel }+e_j)e_j+\mathrm{\parallel }{E}_j\mathrm{\parallel }\mathrm{.}$ The term $E_j$ records the entire defect in the computed squaring.

Let $A\in {\mathbb{C}}^{n\times n}$ with $n\ge 1$, let $\mathrm{\Omega }\subseteq \mathbb{C}$ be open with $\sigma (A)\subseteq \mathrm{\Omega }$, where $\sigma (A)$ is the eigenvalue set from Eigenvalues, eigenvectors, eigenspaces $E_{\lambda }(T)=\ker (T-\lambda I)$, and the spectrum ${\sigma }_F(T)$ of an endomorphism, and let $f:\mathrm{\Omega }\to \mathbb{C}$ be holomorphic in the sense of Complex differentiability at a point, the complex derivative, holomorphic functions, and entire functions. Choose a Jordan representation in the convention of Jordan bases and Jordan canonical forms over the base field, which exists by Jordan form over the base field exists exactly when the characteristic polynomial splits and Fundamental theorem of algebra: every nonconstant complex polynomial has a complex root: $A=Z\mathrm{diag}(J_{d_1}({\lambda }_1),\dots ,J_{d_r}({\lambda }_r))Z^{-1},$ where $J_d(\lambda )=\lambda I_d+N_d$ and $N_d$ has ones on its first superdiagonal and is nilpotent of order at most $d$. Define the Jordan expression $f_J(A):=Z\mathrm{diag}({\sum }_{k=0}^{d_1-1}\frac{f^{(k)}({\lambda }_1)}{k!}{N}_{d_1}^k,\dots ,{\sum }_{k=0}^{d_r-1}\frac{f^{(k)}({\lambda }_r)}{k!}{N}_{d_r}^k)Z^{-1}\mathrm{.}$ The analytic primary matrix function $f(A)$ is this expression, with the same single scalar branch $f$ and its derivatives used in every block at a repeated eigenvalue. Its independence of the chosen Jordan representation is discharged by Jordan, contour, similarity, and normal-functional-calculus formulas agree ↗.

For $m\ge 0$ and $\rho \ge 0$, put $d_m(\rho ):=1-{\sum }_{j=1}^m{a}_{m,j}{\rho }^j,$ with the empty sum equal to zero and $a_{m,j}$ as in Diagonal Padé approximant to the exponential. If $d_m(\rho )>0$, then $q_m$ has no zero in $\mathrm{\mid }z\mathrm{\mid }\le \rho$ and, for every such $z$, $\mid e^z-r_m(z)\mid \le \frac{{\rho }^{2m+1}{e}^{\rho }(m!{)}^2}{(2m)!(2m+1)!d_m(\rho )}\mathrm{.}$ This is a finite-order bound; it makes no adaptive choice of $m$ or scaling.

Let $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, have $\sigma (A)\cap i\mathbb{R}=\mathrm{\varnothing }$, where $\sigma (A)$ is the spectrum from Eigenvalues, eigenvectors, eigenspaces $E_{\lambda }(T)=\ker (T-\lambda I)$, and the spectrum ${\sigma }_F(T)$ of an endomorphism. The scalar function $s(z)=1$ for $\mathrm{Re}z>0$ and $s(z)=-1$ for $\mathrm{Re}z<0$ is locally constant and holomorphic on $\mathbb{C}\setminus i\mathbb{R}$. Define the primary matrix sign by $\mathrm{sign}(A):=s(A)$ using Analytic primary matrix functions by Jordan blocks. This definition fixes the same half-plane branch on every Jordan block belonging to an eigenvalue.

Let $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, and let $f$ be holomorphic on an open set $\mathrm{\Omega }$ containing $\sigma (A)$. The Jordan expression in Analytic primary matrix functions by Jordan blocks is independent of the chosen Jordan representation and satisfies $f(XA{X}^{-1})=Xf(A)X^{-1}$ for every invertible $X$. It equals $p(A)$ for the unique Hermite polynomial $p$ matching the derivatives of $f$ at the Jordan eigenvalue data; in particular, it commutes with $A$. For holomorphic $f,g$ on $\mathrm{\Omega }$, it also satisfies $(fg)(A)=f(A)g(A)$. If $A$ is normal, it equals the spectral functional calculus The spectral functional calculus f(T) for a normal endomorphism.

If $T$ is block upper triangular with diagonal blocks $T_{ii}$ and $f$ is holomorphic on a neighborhood of $\sigma (T)$, then $f(T)$ is block upper triangular and its diagonal blocks are $f(T_{ii})$.

For the contour formula, let $\mathrm{\Gamma }$ be a finite complex cycle with trace in $\mathrm{\Omega }$, null-homologous in $\mathrm{\Omega }$, disjoint from $\sigma (A)$, and with index $n(\mathrm{\Gamma },\lambda )=1$ for every $\lambda \in \sigma (A)$. Then, with the matrix line integral taken entry by entry, $f(A)=\frac{1}{2\pi i}{\int }_{\mathrm{\Gamma }}f(z)(zI-A{)}^{-1}dz\mathrm{.}$

Input a complex matrix $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, an open $\mathrm{\Omega }\subseteq \mathbb{C}$ containing $\sigma (A)$ and a holomorphic function $f:\mathrm{\Omega }\to \mathbb{C}$ specified by exact derivative data of that same $f$ on the spectrum, and an exact supplied Jordan certificate $A=ZJ{Z}^{-1}$ with $J=\mathrm{diag}(J_{d_1}({\lambda }_1),\dots ,J_{d_r}({\lambda }_r))$.

1. Verify the dimensions, exact equality $A=ZJ{Z}^{-1}$, invertibility of $Z$, and the stated Jordan block pattern. If any check fails, return invalid-Jordan-certificate.
2. For each block, form the finite sum $F_i={\sum }_{k=0}^{d_i-1}{f}^{(k)}({\lambda }_i)N_{d_i}^k\mathrm{/}k!$ using the supplied exact derivatives. If a required value is absent, return missing-jet.
3. Return $Z\mathrm{diag}(F_1,\dots ,F_r)Z^{-1}$ with status exact.

All loops and matrix products have finite size. This is an exact-data evaluator; it makes no floating-point stability or complexity claim.

If Finite analytic evaluation from a supplied exact Jordan decomposition returns status exact for the supplied $A$, $f$, and Jordan certificate, then its output is $f(A)$ as defined by Analytic primary matrix functions by Jordan blocks.

Let $S=\mathbb{C}\setminus (-\mathrm{\infty },0]$, and let $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, satisfy $\sigma (A)\subset S$. The analytic primary matrices $L(A)=(\mathrm{Log})(A)$ and $R(A)=(z_{\mathrm{Log}}^{})(A)$ exist and satisfy $\exp (L(A))=A,R(A{)}^2=A\mathrm{.}$ Moreover, $\sigma (L(A))$ lies in the strip $\{w:-\pi <\mathrm{Im}w<\pi \}$, and $\sigma (R(A))$ lies in the open right half-plane. These are the unique matrices $Y$ and $Q$ satisfying respectively $\exp (Y)=A$ and $Q^2=A$ with the stated spectral restrictions; this specifies what “logarithm” and “square root” mean in the uniqueness clause.

Let $\mathrm{\Omega }\subseteq \mathbb{C}$ be open, $f:\mathrm{\Omega }\to \mathbb{C}$ holomorphic, $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, with $\sigma (A)\subset \mathrm{\Omega }$, and $E\in {\mathbb{C}}^{n\times n}$. The map $A\mapsto f(A)$ is real Fréchet differentiable at $A$, and $Df(A)[E]=\frac{1}{2\pi i}{\int }_{\mathrm{\Gamma }}f(z)(zI-A{)}^{-1}E(zI-A{)}^{-1}dz$ for any admissible Cauchy cycle $\mathrm{\Gamma }$ from Jordan, contour, similarity, and normal-functional-calculus formulas agree. Consequently $f\text{\hspace{0.17em}⁣}\left(\left(\begin{array}{cc}A& E\\ 0& A\end{array}\right)\right)=\left(\begin{array}{cc}f(A)& Df(A)[E]\\ 0& f(A)\end{array}\right)\mathrm{.}$ For $f(z)=e^z$, the same derivative also has the formula $D\exp (A)[E]={\int }_0^1{e}^{(1-t)A}E{e}^{tA}dt\mathrm{.}$

Let $A\in {\mathbb{C}}^{n\times n}$ with $n\ge 1$, let $f$ be holomorphic on an open $\mathrm{\Omega }\subseteq \mathbb{C}$ containing $\sigma (A)$, and let $r$ be rational with no pole on $\mathrm{\Omega }$. Let $\mathrm{\Gamma }$ be an admissible Cauchy cycle for $A$ contained in $\mathrm{\Omega }$. Write its finite chain representation as $\mathrm{\Gamma }={\sum }_{\nu }{m}_{\nu }{\gamma }_{\nu }$ with integer coefficients and define its coefficient-weighted length by $L_{\mathrm{\Gamma }}:={\sum }_{\nu }\mathrm{\mid }{m}_{\nu }\mathrm{\mid }L({\gamma }_{\nu })$. The integral and this length use the same signed-chain representation. Suppose ${\sup }_{z\in \mathrm{\Gamma }}\mathrm{\mid }f(z)-r(z)\mathrm{\mid }\le \epsilon ,{\sup }_{z\in \mathrm{\Gamma }}\mathrm{\parallel }(zI-A{)}^{-1}{\mathrm{\parallel }}_{\mathrm{\infty },\mathbb{C}}\le M\mathrm{.}$ Then $\mathrm{\parallel }f(A)-r(A){\mathrm{\parallel }}_{\mathrm{\infty },\mathbb{C}}\le \frac{L_{\mathrm{\Gamma }}}{2\pi }\epsilon M\mathrm{.}$ The resolvent hypothesis may be supplied by the diskwise certificate of Diskwise Neumann certificates bound a complex matrix resolvent. A scalar disk error alone does not give this matrix bound for a nonnormal matrix.

Input exact complex matrices $M\in {\mathbb{C}}^{n\times n}$, $C\in {\mathbb{C}}^{n\times r}$, $n,r\ge 1$, and a nonnegative integer candidate-operation budget $K$. Retain their exact real/imaginary values for the final audit. Rational complex inputs $M,C\in \mathbb{Q}(i)$ are included.

The verification/controller model is an exact real-RAM: its scalar field operations $+,-,\times ,\mathrm{/}$, zero/order comparisons, finite indices and copies act on the exact input/stored real values. Complex controller operations mean their real-component field formulas. The final audit uses the fixed row-major loops for $M\widehat{X}-C$, with inner products summed in increasing index order, and tests all real/imaginary components for zero. These are prescribed finite expression trees/loops, not an unprovided certification oracle. No implementation or bit-cost claim for arbitrary real encodings is made. In a rational adapter, exact rational inputs and a supplied finite format codebook give every stored real code its exact rational value; decoding is finite table lookup. When all required codes/values are supplied, the same audit is finite exact rational arithmetic (binary formats give dyadic values). Missing/invalid adapter data return representation-inconclusive; they are not silently assumed. The default exact real-RAM route requires no rational/dyadic decoder.

Candidate operations use Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Charge one budget unit BEFORE each rounded real primitive, including each component conversion and every primitive inside a called norm, reflector, complex product, subtraction or division. On exhaustion return budget-exhausted with any available partial candidate. Successful exact controller operations, index loops, comparisons, conjugations and copies are separate from this rounded-operation budget; their loops are finite and no total cubic or bit-cost assertion is made here. Each excluded rounded intermediate returns range-failure with its index and partial values. Required stored scalar constants must be admissible, or the candidate route returns candidate-failure.

Instrumented kernel execution. The calls in steps2–3 are macros, not opaque subroutine calls: expand the actual finite loops and scalar expressions of the two named Householder definitions, including the scaled real-component norm, phase-modulus norm, stored-vector squared norm, all sequential real dot operations, real-scalar component multiplication, six-real-operation complex products, additions/subtractions, divisions and square roots. Replace each requested rounded real primitive or component conversion by the following dispatcher. With counter $b$ initially zero, first test $b=K$; if equal, return budget-exhausted immediately without evaluating that primitive. Otherwise increase $b$ by one and execute exactly that one primitive in the source arithmetic model; propagate its first failure. No macro, norm or reflector is executed before its primitives are dispatched, and no budget count is guessed after a completed call. Keep all source loop indices, zero branches, target overwrites, success/failure labels and operation order; a zero-column marker represents the identity and skips reflector denominators and RHS arithmetic. Exact control comparisons/copies remain separate. This is a budgeted execution of the same kernels, not a replacement numerical method. A fully funded execution is identical to the source operation trace. A partial candidate includes an explicit completed-entry mask and is never presented as a complete solution matrix. The exact audit in step4 is reached only after every one of the $nr$ solution entries is successfully computed.

1. Convert $M$ completely before $C$, visiting columns increasingly, rows increasingly within each column, and each real component before its imaginary component, to form stored candidates $R,D$. A failed conversion returns range-failure (or budget-exhausted before an unaffordable conversion).
2. For $k=1,\dots ,n$, call Scaled complex Householder reflector step with explicit range status with target $R_{k:n,k}$ and trailing block $R_{k:n,k+1:n}$. Retain its target overwrite and trailing update, and store $v,\widehat{\tau }$ and the represented factor. Propagate breakdown, range-failure and budget exhaustion. A zero-column identity branch copies the blocks. Replay this SAME stored reflector on the active RHS panel $D_{k:n,:}$ in formation order, without any target overwrite, using the separate-block kernel of Scaled real and complex Householder and Givens update arithmetic: for each RHS column $j$ compute $d_j=\mathrm{fl}\text{\hspace{0.17em}⁣}\left({\sum }_{i=k}^n{\bar{v}}_i{D}_{ij}\right),p_i=\mathrm{fl}(\widehat{\tau }{v}_i),z_{ij}=\mathrm{fl}(p_i{d}_j),D_{ij}\leftarrow \mathrm{fl}(D_{ij}-z_{ij})\mathrm{.}$ Use the source's increasing indices, first-term-copy sequential dots, componentwise real-coefficient multiplication and declared complex products. Finish all $d_j$, then recompute all $p_i$ from the stored $v,\widehat{\tau }$, then all $z_{ij}$, then all subtractions, in the source's column-major order. No product from a previous block call is reused. Identity reflectors copy $D$. Do not feed nonstored normalized reflector entries to floating operations. After the column step, if $R_{kk}=0$, return zero-pivot.
3. For each RHS $j$ in increasing order and each row $i=n,\dots ,1$, start $s=D_{ij}$ and for $\mathrm{\ell }=i+1,\dots ,n$ set $s\leftarrow \mathrm{fl}(s-\mathrm{fl}(R_{i\mathrm{\ell }}{\widehat{X}}_{\mathrm{\ell }j}))$. This is the rounded backward recurrence of Forward and backward substitution for triangular linear systems. To divide $s=a+ib$ by the nonzero stored pivot $R_{ii}=c+id$, use ONLY the following rounded real primitives in the displayed order: $t=\mathrm{fl}(\mathrm{fl}(cc)+\mathrm{fl}(dd)),p=\mathrm{fl}(\mathrm{fl}(ac)+\mathrm{fl}(bd)),q=\mathrm{fl}(\mathrm{fl}(bc)-\mathrm{fl}(ad)),$ ${\widehat{X}}_{ij}=\mathrm{fl}(p\mathrm{/}t)+i\mathrm{fl}(q\mathrm{/}t)\mathrm{.}$ Reject $t\le 0$ as range-failure; every other nonzero excluded intermediate also fails under the scalar model. This deliberately conservative routine can fail on range even for a nonzero pivot; it is not a new complex-division primitive. Other invalid candidate data return candidate-failure.
4. Compute $M\widehat{X}-C$ in the declared exact controller on the original inputs and exact stored values. If it is exactly zero, return verified-solution with $\widehat{X}$; otherwise return candidate-unverified with $\widehat{X}$ and the exact nonzero residual. Use the rational adapter only when its actual encoded-data preconditions hold.