Skip to content

[feat] squareroot matrix function - #278

Draft
lkdvos wants to merge 6 commits into
mainfrom
ld-squareroot
Draft

[feat] squareroot matrix function#278
lkdvos wants to merge 6 commits into
mainfrom
ld-squareroot

Conversation

@lkdvos

@lkdvos lkdvos commented Aug 26, 2026

Copy link
Copy Markdown
Member

This is an attempt to move forward with #261 in smaller pieces, mostly since I got stuck a bit with the various different domain tolerances. Additionally I actually wanted to try out the actual square root dedicated algorithms to have something to go beyond the normal linear algebra routines, so here is that progress. What's nice is that this approach actually already works for non-blas paths, and seems to again outperform the linearalgebra routines.

To do

  • clean up PR
  • actually perform benchmarks
  • review changes properly

lkdvos and others added 6 commits August 15, 2026 05:15
The exponential tests lived in three places: `common/exponential.jl` for the
BLAS floats and one file each in the `genericschur` and `genericlinearalgebra`
groups for the generic ones. Move the assertions into the shared TestSuite and
leave a single driver that covers every element type, which removes both
generic groups.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`squareroot` computes the principal square root, through `LinearAlgebra`, an
(hermitian) eigenvalue decomposition, or the `Diagonal` fast path. The scalar
type of the output matches the input, so a real matrix whose spectrum reaches
the negative real axis raises a `DomainError`; eigenvalues that are negative
only within `domain_atol` are clamped onto zero instead.

`MatrixFunctionViaEig` and `MatrixFunctionViaEigh` gain a `domain_atol` field
for this, and `src/implementations/matrixfunctions.jl` collects the input
handling, reconstruction and domain helpers that the other matrix functions
will share.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The assertions are invariants — `sqrt(A)^2 ≈ A`, hermiticity, the elementwise
action on the spectrum — so that they also apply on GPU and to downstream array
types, with a host-only cross-check against `LinearAlgebra`. The domain tests
prescribe the full spectrum, and probe a genuinely negative eigenvalue, one that
is negative by roundoff, and one in between that only an explicit `domain_atol`
admits.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Describe the shared matrix function algorithms once, and add a section on the
domain of `squareroot`: what `domain_atol` clamps, why an accepted result can be
off by more than the tolerance, and why the tolerance of `MatrixFunctionViaLA`
measures something else entirely.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant