Skip to content

Commit

Permalink
Fix log_quasitriu for internal scaling s=0 (JuliaLang#56311)
Browse files Browse the repository at this point in the history
This PR is a potential fix for #54833.

## Description
The function
https://github.com/JuliaLang/julia/blob/2a06376c18afd7ec875335070743dcebcd85dee7/stdlib/LinearAlgebra/src/triangular.jl#L2220
computes $\boldsymbol{A}^{\dfrac{1}{2^s}} - \boldsymbol{I}$ for a
real-valued $2\times 2$ matrix $\boldsymbol{A}$ using Algorithm 5.1 in
[R1]. However, the algorithm in [R1] as well as the above function do
not handle the case $s=0.$ This fix extends the function to compute
$\boldsymbol{A}^{\dfrac{1}{2^s}} - \boldsymbol{I} \Bigg|_{s=0} =
\boldsymbol{A} - \boldsymbol{I}.$

## Checklist
- [X] Fix code: `stdlib\LinearAlgebra\src\triangular.jl` in function
`_sqrt_pow_diag_block_2x2!(A, A0, s)`.
- [X] Add test case: `stdlib\LinearAlgebra\test\triangular.jl`.
- [X] Update `NEWS.md`.
- [X] Testing and self review.

|  Tag  | Reference |
| --- | --- |
| <nobr>[R1]</nobr> | Al-Mohy, Awad H. and Higham, Nicholas J. "Improved
Inverse Scaling and Squaring Algorithms for the Matrix Logarithm", 2011,
url: https://eprints.maths.manchester.ac.uk/1687/1/paper11.pdf |

---------

Co-authored-by: Daniel Karrasch <[email protected]>
Co-authored-by: Oscar Smith <[email protected]>
  • Loading branch information
3 people authored Oct 28, 2024
1 parent 87ecf8f commit 2cdfe06
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 0 deletions.
5 changes: 5 additions & 0 deletions stdlib/LinearAlgebra/src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2218,6 +2218,11 @@ end
# SIAM J. Sci. Comput., 34(4), (2012) C153–C169. doi: 10.1137/110852553
# Algorithm 5.1
Base.@propagate_inbounds function _sqrt_pow_diag_block_2x2!(A, A0, s)
if iszero(s)
A[1,1] -= 1
A[2,2] -= 1
return A
end
_sqrt_real_2x2!(A, A0)
if isone(s)
A[1,1] -= 1
Expand Down
10 changes: 10 additions & 0 deletions stdlib/LinearAlgebra/test/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1386,4 +1386,14 @@ end
end
end

@testset "log_quasitriu with internal scaling s=0 (issue #54833)" begin
M = [0.9949357359852791 -0.015567763143324862 -0.09091193493947397 -0.03994428739762443 0.07338356301650806;
0.011813655598647289 0.9968988574699793 -0.06204555000202496 0.04694097614450692 0.09028834462782365;
0.092737943594701 0.059546719185135925 0.9935850721633324 0.025348893985651405 -0.018530261590167685;
0.0369187299165628 -0.04903571106913449 -0.025962938675946543 0.9977767446862031 0.12901494726320517;
0.0 0.0 0.0 0.0 1.0]

@test exp(log(M)) M
end

end # module TestTriangular

0 comments on commit 2cdfe06

Please sign in to comment.