feat: `squareroot`, `logarithm` and `power` by lkdvos · Pull Request #261 · QuantumKitHub/MatrixAlgebraKit.jl

lkdvos · 2026-07-27T22:58:56Z

Rebased onto main, which now contains #262 (the Pkg-workspace test migration). This PR adds the test/matrixfunctions group that #262's description deferred here, and updates test/README.md accordingly.

Following up on the exponential implementation, with the goal of getting QuantumKitHub/TensorKit.jl#493 to work.

I do have the intention of adding better algorithms for the squareroot and logarithm implementations, but I figured it might be nice to avoid doing that in one huge PR and to keep things moving by defining the interface first.

Still to do:

Refactor tests into testsuite
Decide if we want to restructure the folders to separate decompositions and matrix functions
See if I can refactor some of these implementations to share a bit more code
Add tests for GPU support

Notes for review

A few things worth knowing, mostly about behaviour that was discovered rather than designed:

power(A, 0) and power(A, 1) short-circuit to I and A before any decomposition, so they
now apply to any square matrix whatever the algorithm — including input the algorithm would
otherwise reject, e.g. a non-hermitian matrix with MatrixFunctionViaEigh. A^0 = I and
A^1 = A hold regardless of spectrum, and Base.^ makes no hermiticity demand either.
An integral Float exponent is still recognised at runtime, via isinteger(p) rather than
dispatch on typeof(p). This matches Base, whose ^(::AbstractMatrix, ::Real) does
isinteger(p) && return integerpow(A, p); dispatching on the type instead would make
power(A, 3.0) throw a DomainError for a matrix with a negative eigenvalue while A^3.0
succeeds.
exponential of a complex hermitian matrix is now exactly hermitian, which it was not before:
the real branch built a symmetric product but the complex branch used a plain V exp(D) V' with
no closing projection. A complex τ is still left unprojected, because exp(τA) is hermitian
only for real τ.
The domain rules genuinely differ per function, and the tests encode that rather than
smoothing it over: a roundoff-negative eigenvalue clamps for squareroot and positive fractional
power, but still throws for logarithm, where clamping onto zero leaves the matrix singular.
MatrixFunctionViaLA inspects only the realness of its result, so it accepts singular input that
the spectrum-based algorithms reject.

Two pre-existing issues on main were found while testing and are deliberately not addressed
here:

MatrixFunctionViaEig is not GPU-clean — its lu!-based solve scalar-indexes — and
MatrixFunctionViaLA calls LAPACK on device memory. Device support for squareroot, logarithm
and power is therefore hermitian-via-eigh plus the Diagonal fast path; exponential escapes
this through the native MatrixFunctionViaTaylor.
power(A, -1) on a numerically singular matrix returns garbage rather than throwing for the
eig-based algorithms, because any(iszero, λ) is an exact test against computed eigenvalues.
The Diagonal path throws SingularException correctly.

codecov · 2026-07-28T02:04:44Z

Codecov Report

❌ Patch coverage is 84.54106% with 32 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/interface/logarithm.jl	0.00%	7 Missing ⚠️
src/interface/power.jl	0.00%	7 Missing ⚠️
src/interface/squareroot.jl	0.00%	7 Missing ⚠️
ext/MatrixAlgebraKitGenericSchurExt.jl	0.00%	3 Missing ⚠️
src/implementations/logarithm.jl	94.28%	2 Missing ⚠️
src/implementations/power.jl	95.83%	2 Missing ⚠️
src/implementations/squareroot.jl	94.11%	2 Missing ⚠️
src/implementations/exponential.jl	80.00%	1 Missing ⚠️
src/implementations/matrixfunctions.jl	97.87%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/MatrixAlgebraKit.jl	`100.00% <ø> (ø)`
src/common/defaults.jl	`66.66% <100.00%> (+2.38%)`	⬆️
src/implementations/polar.jl	`96.05% <ø> (-0.18%)`	⬇️
src/implementations/projections.jl	`92.63% <100.00%> (+0.58%)`	⬆️
src/interface/matrixfunctions.jl	`93.75% <100.00%> (-6.25%)`	⬇️
src/implementations/exponential.jl	`96.04% <80.00%> (-1.13%)`	⬇️
src/implementations/matrixfunctions.jl	`97.87% <97.87%> (ø)`
src/implementations/logarithm.jl	`94.28% <94.28%> (ø)`
src/implementations/power.jl	`95.83% <95.83%> (ø)`
src/implementations/squareroot.jl	`94.11% <94.11%> (ø)`
... and 4 more

... and 4 files with indirect coverage changes

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Adds the PR link to the existing entries, which every other entry in the file carries, plus the two user-visible changes made since: the trivial exponents of `power` short-circuiting the decomposition, and `exponential` of a complex hermitian matrix now being exactly hermitian. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Jutho · 2026-07-29T09:22:14Z

+    # so a complex result with a real output signals a genuine domain violation
+    logAc = LinearAlgebra.log(A)
+    if eltype(logAc) <: Complex && !(eltype(logA) <: Complex)
+        throw(_realness_domainerror(logarithm!))


Do we need a tolerance check on the imaginary parts as for the case of power, or is the situation different here?

Jutho · 2026-07-29T09:34:37Z

+struct MatrixFunctionViaEigh{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm
    eigh_alg::A
+    domain_atol::T
+end


I guess changing this struct definition makes this a breaking change?

I guess you are right... There might be an argument to be made for hardcoding the tolerance to Float64, not adding the type parameter, and then this looks non-breaking to me?

In the following round of breaking changes, we could also consider making this actual @algdefs, which would have solved this problem...

Jutho · 2026-07-29T09:44:35Z

I did not go through all the test changes, also because this contains stuff from the other PR. The changes to src, ext and docs look good to me, though for the latter, it would be useful if we can also have a github action that builds a docs preview and put the links here, like in some of the other repos.

Jutho · 2026-07-29T09:51:49Z

One thing I just noted in the docs is that we have a section on "Matrix Properties" and then one on "Matrix functions", i.e. we have a capitalization inconsistency 😄 .

Following the `exponential` template, each function supports `MatrixFunctionViaLA` (wrapping the LinearAlgebra implementation), `MatrixFunctionViaEig`/`MatrixFunctionViaEigh` (eigendecomposition-based, also serving generic eltypes via GenericSchur/GenericLinearAlgebra), and `DiagonalAlgorithm`. `power(A, p)` uses the two-argument `@functiondef` form with an exponent-type-based split between integer and fractional powers. The output scalar type always matches the input: out-of-domain eigenvalues (negative real axis for real input, zero eigenvalues for `logarithm` and negative fractional powers) throw a `DomainError`, while eigenvalues that violate the domain within a tolerance `domain_atol` (new keyword of `MatrixFunctionViaEig(h)`, defaulting to the new `default_domain_atol`) are clamped as rounding artifacts.

Mirrors the `exponential` test layout: BLAS floats against the LinearAlgebra reference, hermitian fast paths, domain-error and tolerance-clamping edge cases, `Diagonal` inputs over generic float types, and BigFloat coverage via the GenericSchur and GenericLinearAlgebra test directories.

Move `_mul_herm!` from `implementations/polar.jl` to `implementations/projections.jl` and use it for the symmetric-product kernels of `squareroot` and fractional `power` via `MatrixFunctionViaEigh`, so hermitian output is guaranteed through the same mechanism as the polar decomposition (BLAS `herk`/`syrk`, explicit projection otherwise, GPU overrides). Rewrite the `MatrixFunctionViaLA` kernels in the style of the `exponential` implementation: an explicit eltype branch with a `real.(...)` assignment instead of the bespoke `_copy_result!` helper. `squareroot`/`logarithm` use a strict complex-result check (LinearAlgebra returns real results whenever the principal value is real), while `power` tolerates rounding-level imaginary components since fractional powers are computed in complex arithmetic upstream.

…nels In the `MatrixFunctionViaEig(h)` implementations of `squareroot`, `logarithm` and `power`, apply the scalar function to the eigenvalue matrix by calling the function's own `DiagonalAlgorithm` kernel (which already contains the domain checks and tolerance clamping), mirroring how `exponential!` recurses into `exponential!((τ, D), D)`. The eig kernels only retain the negative-real-axis scan of the complex eigenvalues for real input, which has no `Diagonal` counterpart.

Add a `_clamp_domain_eigenvalues!(D::Diagonal, domain_atol)` convenience method that derives the default tolerance itself, so the real-input branch of the `MatrixFunctionViaEig` kernels reduces to a one-line preamble before delegating to the `DiagonalAlgorithm` kernel.

Replace the scalar-indexing loops in `_clamp_domain_eigenvalues!` and `_check_nonzero_eigenvalues` with reductions and broadcasts, so the helpers work on GPU arrays. The `DomainError` now reports the worst offending eigenvalue instead of the first one encountered.

The `lazy"..."` literal introduced for the negative-eigenvalue error does not process the `\` line continuation, so the rendered message contained a literal backslash, newline and indentation. Use ordinary string concatenation, as elsewhere in this file, and report the tolerance under its user-facing name `domain_atol` rather than `atol`. Restore the distinction between the real ("a negative real eigenvalue") and complex ("an eigenvalue on the negative real axis") variants by passing the description to the shared helper, and give `_check_nonzero_eigenvalues` the same `@noinline` throw helper so both domain checks stay free of error-path code on the GPU. Thread `domain_atol` back into the diagonal kernels invoked by the `squareroot` eigendecomposition paths, and let the `DiagonalAlgorithm` kernel reuse `_clamp_domain_eigenvalues!` again so it matches `power!`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

The domain helpers were rewritten as reductions and broadcasts so they would run on GPU arrays, but nothing exercised that. Add `test_{squareroot,logarithm,power}_gpu`, modelled on the existing `test_exponential_gpu`, comparing device results against the CPU reference for the `DiagonalAlgorithm` fast path and the `MatrixFunctionViaEigh` path. Both domain outcomes are covered, since the clamp and the throw are precisely what the scalar loops were replaced by: roundoff-scale negative eigenvalues must clamp, genuinely negative ones must throw a `DomainError`, and for `logarithm` and negative fractional `power` a zero eigenvalue must throw as well. Scalar indexing anywhere along these paths errors under the GPUArrays guard, so a regression to the previous loops would fail these tests. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

`map_diagonal(f, src)` is exactly `diagonal(f.(diagview(src)))`, so the three `copy_input` methods can use it instead of constructing a `Diagonal` directly. This keeps them consistent with the rest of the package, which routes diagonal construction through `diagonal` rather than the `LinearAlgebra` type. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

The four matrix functions were tested by hand-written, CPU-only top-level files plus per-backend duplicates under `test/genericschur/` and `test/genericlinearalgebra/`, while every decomposition lives in `test/testsuite/` as array-type-parameterized functions driven by a thin `test/<area>/<name>.jl`. Close that gap: `test/testsuite/matrixfunctions/` now holds the bodies and `test/matrixfunctions/` the drivers. Assertions become invariants rather than comparisons against `LinearAlgebra`, so one body serves host, CUDA, ROCm and downstream array types: `sqrtA * sqrtA ≈ A`, `exponential(logA) ≈ A`, `exp(A) * exp(-A) ≈ I`, and elementwise spectrum mapping via `eigh_vals`. Reference cross-checks against `LinearAlgebra` are kept, but confined to `test_*_reference`, which the drivers only call for host arrays. GPU coverage now falls out of the array-type parameterization, so the bespoke `test_*_gpu` functions are gone. Three new instantiators (`instantiate_offaxis_matrix`, `instantiate_posdef_matrix`, `instantiate_hermitian_spectrum`) give the controlled spectra these functions need, reusing `instantiate_unitary` and the `Diagonal` handling already in the suite. The domain rules differ per function and are now encoded rather than approximated: a roundoff-negative eigenvalue clamps for `squareroot` and positive fractional `power`, but still throws for `logarithm`, where clamping onto zero leaves the matrix singular. `MatrixFunctionViaLA` inspects only the realness of its result, so it accepts singular input that the spectrum-based algorithms reject; the `eigh`-based algorithms reject a negative eigenvalue whatever the scalar type, since the root of a hermitian matrix with a negative eigenvalue is not hermitian. Also drops the deprecated `LAPACK_Simple`/`LAPACK_QRIteration`/ `GS_QRIteration`/`GLA_QRIteration` spellings, which accounted for the deprecation warnings these tests emitted. Generic eltypes name their driver explicitly (`QRIteration(; driver = GS())`), since with both generic packages loaded `default_driver` resolves `QRIteration` to GLA, which provides no `geev!`. Full suite: 93826 passing, up from 90917. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

…tions `squareroot`, `logarithm` and `power` all map a square matrix to a matrix of the same size and scalar type, so their `copy_input`, `check_input` and `initialize_output` contracts were identical modulo the function name -- six methods repeated three times over. Fold the bodies into `_matrixfunction_copy_input` and `_matrixfunction_check_input`, whose own dispatch absorbs the `Diagonal` and `DiagonalAlgorithm` specializations, so each function keeps just the one line that names it. The three `MatrixFunctionViaEig` kernels shared an identical skeleton, differing only in which scalar function is applied to the eigenvalues: solve `V f(D) V⁻¹` and, for real input, discard the imaginary part left by the complex decomposition. That becomes `_apply_eig!`, taking the already transformed eigenvalues. Similarly `_apply_eigh!` captures the `V f(D) V'`-then-project shape shared by `logarithm` and integer `power`, and pairs with the existing `_mul_herm!` for the cases where a symmetric product makes the result hermitian by construction instead. Deliberately left alone: the `MatrixFunctionViaLA` kernels, whose bodies look alike but must not be merged -- `_copy_result!` did exactly that and was removed in 2698e26, because `squareroot`/`logarithm` need a strict complex-result check while `power` has to tolerate rounding-level imaginary components. `exponential` also stays out, as its real/complex branch keys on the scalar `τ` as well as the matrix eltype. Behaviour is unchanged: the matrix-function tests report the same 4701 passing as before the refactor. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

`MatrixFunctionViaEigh` built the real case as the symmetric product `VexpD * transpose(VexpD)`, exact by construction, but the complex case as a plain `VexpD * V'` with no closing projection, so the result was only approximately hermitian. `squareroot`, `logarithm` and `power` are all exact for both scalar types, leaving `exponential` the odd one out. Since `eigh_full!` returns real eigenvalues, `exp(τD/2)` is self-adjoint whenever `τ` is real, and `exp(τA) = (V exp(τD/2)) * (V exp(τD/2))'` is then hermitian by construction. Route the complex case through `_mul_herm!`, the same mechanism the other three use. This needed a second change to take effect for `exponential(A)` itself: the implicit scalar was `one(eltype(A))`, i.e. `1.0 + 0.0im` for a complex matrix, which sent a plain `exp(A)` down the complex-`τ` branch even though the scalar was real. Use `one(real(eltype(A)))` instead. The same edit for `MatrixFunctionViaEig` is behaviour-neutral -- its branch condition already excluded a complex `A` -- but avoids needless complex-scalar arithmetic. A genuinely complex `τ` is still left unprojected, which is not an oversight: `exp(τA)` of a hermitian `A` is hermitian only for real `τ`, and projecting would silently return a different matrix. Verified against `LinearAlgebra.exp(Hermitian(A))` for `Float32`, `Float64`, `ComplexF32` and `ComplexF64` with both `QRIteration` and `DivideAndConquer`, and for `Diagonal` with complex storage. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Adds the PR link to the existing entries, which every other entry in the file carries, plus the two user-visible changes made since: the trivial exponents of `power` short-circuiting the decomposition, and `exponential` of a complex hermitian matrix now being exactly hermitian. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

The README added in #262 describes the layout this branch changes: the matrix exponential has moved out of `common`, the `genericschur` and `genericlinearalgebra` groups are gone now that each driver covers the generic element types alongside the BLAS floats, and `matrixfunctions` is a new group. Also records why the `eig`-based algorithms have to name their driver explicitly for those element types, and lists the spectrum-controlled instantiators the matrix functions need. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

`pkgdocs.julialang.org/dev/workspaces/` is a 404 -- workspaces have no page of their own, they are a section of the `Project.toml` reference. Point at `v1/toml-files/#Workspaces` instead: the anchor resolves to "The `[workspace]` section", and the versioned channel will not move the way the `dev` docs did. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>