Skip to content

Commit

Permalink
Don't try cholesky for Hermitian in factorize (#54775)
Browse files Browse the repository at this point in the history
  • Loading branch information
dkarrasch authored Jun 13, 2024
1 parent c8a90cc commit 1505b07
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 21 deletions.
26 changes: 9 additions & 17 deletions stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1387,15 +1387,14 @@ end
factorize(A)
Compute a convenient factorization of `A`, based upon the type of the input matrix.
`factorize` checks `A` to see if it is symmetric/triangular/etc. if `A` is passed
as a generic matrix. `factorize` checks every element of `A` to verify/rule out
each property. It will short-circuit as soon as it can rule out symmetry/triangular
structure. The return value can be reused for efficient solving of multiple
systems. For example: `A=factorize(A); x=A\\b; y=A\\C`.
If `A` is passed as a generic matrix, `factorize` checks to see if it is
symmetric/triangular/etc. To this end, `factorize` may check every element of `A` to
verify/rule out each property. It will short-circuit as soon as it can rule out
symmetry/triangular structure. The return value can be reused for efficient solving
of multiple systems. For example: `A=factorize(A); x=A\\b; y=A\\C`.
| Properties of `A` | type of factorization |
|:---------------------------|:-----------------------------------------------|
| Positive-definite | Cholesky (see [`cholesky`](@ref)) |
| Dense Symmetric/Hermitian | Bunch-Kaufman (see [`bunchkaufman`](@ref)) |
| Sparse Symmetric/Hermitian | LDLt (see [`ldlt`](@ref)) |
| Triangular | Triangular |
Expand All @@ -1406,9 +1405,6 @@ systems. For example: `A=factorize(A); x=A\\b; y=A\\C`.
| General square | LU (see [`lu`](@ref)) |
| General non-square | QR (see [`qr`](@ref)) |
If `factorize` is called on a Hermitian positive-definite matrix, for instance, then `factorize`
will return a Cholesky factorization.
# Examples
```jldoctest
julia> A = Array(Bidiagonal(fill(1.0, (5, 5)), :U))
Expand All @@ -1427,8 +1423,9 @@ julia> factorize(A) # factorize will check to see that A is already factorized
⋅ ⋅ ⋅ 1.0 1.0
⋅ ⋅ ⋅ ⋅ 1.0
```
This returns a `5×5 Bidiagonal{Float64}`, which can now be passed to other linear algebra functions
(e.g. eigensolvers) which will use specialized methods for `Bidiagonal` types.
This returns a `5×5 Bidiagonal{Float64}`, which can now be passed to other linear algebra
functions (e.g. eigensolvers) which will use specialized methods for `Bidiagonal` types.
"""
function factorize(A::AbstractMatrix{T}) where T
m, n = size(A)
Expand Down Expand Up @@ -1490,12 +1487,7 @@ function factorize(A::AbstractMatrix{T}) where T
return UpperTriangular(A)
end
if herm
cf = cholesky(A; check = false)
if cf.info == 0
return cf
else
return factorize(Hermitian(A))
end
return factorize(Hermitian(A))
end
if sym
return factorize(Symmetric(A))
Expand Down
7 changes: 3 additions & 4 deletions stdlib/LinearAlgebra/test/cholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -61,16 +61,15 @@ end

# Test of symmetric pos. def. strided matrix
apd = a'*a
@inferred cholesky(apd)
capd = factorize(apd)
capd = @inferred cholesky(apd)
r = capd.U
κ = cond(apd, 1) #condition number

unary_ops_tests(apd, capd, ε*κ*n)
if eltya != Int
@test Factorization{eltya}(capd) === capd
if eltya <: Real
@test Array(Factorization{complex(eltya)}(capd)) Array(factorize(complex(apd)))
@test Array(Factorization{complex(eltya)}(capd)) Array(cholesky(complex(apd)))
@test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya)
end
end
Expand Down Expand Up @@ -358,7 +357,7 @@ end
0.11192755545114 - 0.1603741874112385im 0.8439562576196216 + 1.0850814110398734im
-1.0568488936791578 - 0.06025820467086475im 0.12696236014017806 - 0.09853584666755086im]
cholesky(Hermitian(apd, :L), RowMaximum()) \ b
r = factorize(apd).U
r = cholesky(apd).U
E = abs.(apd - r'*r)
ε = eps(abs(float(one(ComplexF32))))
n = 10
Expand Down

0 comments on commit 1505b07

Please sign in to comment.