Skip to content

Commit

Permalink
Add matrix symmetrizing functions (#31836)
Browse files Browse the repository at this point in the history
Co-authored-by: Steven G. Johnson <[email protected]>
Co-authored-by: Alex Arslan <[email protected]>
Co-authored-by: Daniel Karrasch <[email protected]>
  • Loading branch information
3 people authored Feb 15, 2023
1 parent 113c2f3 commit faf7954
Show file tree
Hide file tree
Showing 5 changed files with 140 additions and 35 deletions.
6 changes: 4 additions & 2 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -43,14 +43,16 @@ Standard library changes

#### Package Manager

- "Package Extensions": support for loading a piece of code based on other
* "Package Extensions": support for loading a piece of code based on other
packages being loaded in the Julia session.
This has similar applications as the Requires.jl package but also
supports precompilation and setting compatibility.
- `Pkg.precompile` now accepts `timing` as a keyword argument which displays per package timing information for precompilation (e.g. `Pkg.precompile(timing=true)`)
* `Pkg.precompile` now accepts `timing` as a keyword argument which displays per package timing information for precompilation (e.g. `Pkg.precompile(timing=true)`)

#### LinearAlgebra

* New functions `hermitianpart` and `hermitianpart!` for extracting the Hermitian
(real symmetric) part of a matrix ([#31836]).

#### Printf
* Format specifiers now support dynamic width and precision, e.g. `%*s` and `%*.*g` ([#40105]).
Expand Down
2 changes: 2 additions & 0 deletions stdlib/LinearAlgebra/docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -467,6 +467,8 @@ Base.copy(::Union{Transpose,Adjoint})
LinearAlgebra.stride1
LinearAlgebra.checksquare
LinearAlgebra.peakflops
LinearAlgebra.hermitianpart
LinearAlgebra.hermitianpart!
```

## Low-level matrix operations
Expand Down
2 changes: 2 additions & 0 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,8 @@ export
eigvecs,
factorize,
givens,
hermitianpart,
hermitianpart!,
hessenberg,
hessenberg!,
isdiag,
Expand Down
131 changes: 98 additions & 33 deletions stdlib/LinearAlgebra/src/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,34 +17,45 @@ end
Construct a `Symmetric` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`)
triangle of the matrix `A`.
`Symmetric` views are mainly useful for real-symmetric matrices, for which
specialized algorithms (e.g. for eigenproblems) are enabled for `Symmetric` types.
More generally, see also [`Hermitian(A)`](@ref) for Hermitian matrices `A == A'`, which
is effectively equivalent to `Symmetric` for real matrices but is also useful for
complex matrices. (Whereas complex `Symmetric` matrices are supported but have few
if any specialized algorithms.)
To compute the symmetric part of a real matrix, or more generally the Hermitian part `(A + A') / 2` of
a real or complex matrix `A`, use [`hermitianpart`](@ref).
# Examples
```jldoctest
julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4]
5×5 Matrix{Int64}:
1 0 2 0 3
0 4 0 5 0
6 0 7 0 8
0 9 0 1 0
2 0 3 0 4
julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
1 2 3
4 5 6
7 8 9
julia> Supper = Symmetric(A)
5×5 Symmetric{Int64, Matrix{Int64}}:
1 0 2 0 3
0 4 0 5 0
2 0 7 0 8
0 5 0 1 0
3 0 8 0 4
3×3 Symmetric{Int64, Matrix{Int64}}:
1 2 3
2 5 6
3 6 9
julia> Slower = Symmetric(A, :L)
5×5 Symmetric{Int64, Matrix{Int64}}:
1 0 6 0 2
0 4 0 9 0
6 0 7 0 3
0 9 0 1 0
2 0 3 0 4
3×3 Symmetric{Int64, Matrix{Int64}}:
1 4 7
4 5 8
7 8 9
julia> hermitianpart(A)
3×3 Hermitian{Float64, Matrix{Float64}}:
1.0 3.0 5.0
3.0 5.0 7.0
5.0 7.0 9.0
```
Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if `A == transpose(A)`).
Note that `Supper` will not be equal to `Slower` unless `A` is itself symmetric (e.g. if
`A == transpose(A)`).
"""
function Symmetric(A::AbstractMatrix, uplo::Symbol=:U)
checksquare(A)
Expand Down Expand Up @@ -99,25 +110,33 @@ end
Construct a `Hermitian` view of the upper (if `uplo = :U`) or lower (if `uplo = :L`)
triangle of the matrix `A`.
To compute the Hermitian part of `A`, use [`hermitianpart`](@ref).
# Examples
```jldoctest
julia> A = [1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4];
julia> A = [1 2+2im 3-3im; 4 5 6-6im; 7 8+8im 9]
3×3 Matrix{Complex{Int64}}:
1+0im 2+2im 3-3im
4+0im 5+0im 6-6im
7+0im 8+8im 9+0im
julia> Hupper = Hermitian(A)
5×5 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
1+0im 0+0im 2+2im 0+0im 3-3im
0+0im 4+0im 0+0im 5+0im 0+0im
2-2im 0+0im 7+0im 0+0im 8+8im
0+0im 5+0im 0+0im 1+0im 0+0im
3+3im 0+0im 8-8im 0+0im 4+0im
3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
1+0im 2+2im 3-3im
2-2im 5+0im 6-6im
3+3im 6+6im 9+0im
julia> Hlower = Hermitian(A, :L)
5×5 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
1+0im 0+0im 6+6im 0+0im 2-2im
0+0im 4+0im 0+0im 9+0im 0+0im
6-6im 0+0im 7+0im 0+0im 3+3im
0+0im 9+0im 0+0im 1+0im 0+0im
2+2im 0+0im 3-3im 0+0im 4+0im
3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}:
1+0im 4+0im 7+0im
4+0im 5+0im 8-8im
7+0im 8+8im 9+0im
julia> hermitianpart(A)
3×3 Hermitian{ComplexF64, Matrix{ComplexF64}}:
1.0+0.0im 3.0+1.0im 5.0-1.5im
3.0-1.0im 5.0+0.0im 7.0-7.0im
5.0+1.5im 7.0+7.0im 9.0+0.0im
```
Note that `Hupper` will not be equal to `Hlower` unless `A` is itself Hermitian (e.g. if `A == adjoint(A)`).
Expand Down Expand Up @@ -863,3 +882,49 @@ for func in (:log, :sqrt)
end
end
end

"""
hermitianpart(A, uplo=:U) -> Hermitian
Return the Hermitian part of the square matrix `A`, defined as `(A + A') / 2`, as a
[`Hermitian`](@ref) matrix. For real matrices `A`, this is also known as the symmetric part
of `A`; it is also sometimes called the "operator real part". The optional argument `uplo` controls the corresponding argument of the
[`Hermitian`](@ref) view. For real matrices, the latter is equivalent to a
[`Symmetric`](@ref) view.
See also [`hermitianpart!`](@ref) for the corresponding in-place operation.
!!! compat "Julia 1.10"
This function requires Julia 1.10 or later.
"""
hermitianpart(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart(A), uplo)

"""
hermitianpart!(A, uplo=:U) -> Hermitian
Overwrite the square matrix `A` in-place with its Hermitian part `(A + A') / 2`, and return
[`Hermitian(A, uplo)`](@ref). For real matrices `A`, this is also known as the symmetric
part of `A`.
See also [`hermitianpart`](@ref) for the corresponding out-of-place operation.
!!! compat "Julia 1.10"
This function requires Julia 1.10 or later.
"""
hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) = Hermitian(_hermitianpart!(A), uplo)

_hermitianpart(A::AbstractMatrix) = _hermitianpart!(copy_similar(A, Base.promote_op(/, eltype(A), Int)))
_hermitianpart(a::Number) = real(a)

function _hermitianpart!(A::AbstractMatrix)
require_one_based_indexing(A)
n = checksquare(A)
@inbounds for j in 1:n
A[j, j] = _hermitianpart(A[j, j])
for i in 1:j-1
A[i, j] = val = (A[i, j] + adjoint(A[j, i])) / 2
A[j, i] = adjoint(val)
end
end
return A
end
34 changes: 34 additions & 0 deletions stdlib/LinearAlgebra/test/symmetric.jl
Original file line number Diff line number Diff line change
Expand Up @@ -790,4 +790,38 @@ end
end
end

@testset "hermitian part" begin
for T in [Float32, Complex{Float32}, Int32, Rational{Int32},
Complex{Int32}, Complex{Rational{Int32}}]
f, f!, t = hermitianpart, hermitianpart!, T <: Real ? transpose : adjoint
X = T[1 2 3; 4 5 6; 7 8 9]
T <: Complex && (X .+= im .* X)
Xc = copy(X)
Y = (X + t(X)) / 2
U = f(X)
L = f(X, :L)
@test U isa Hermitian
@test L isa Hermitian
@test U.uplo == 'U'
@test L.uplo == 'L'
@test U == L == Y
if T <: AbstractFloat || real(T) <: AbstractFloat
HU = f!(X)
@test HU == Y
@test triu(X) == triu(Y)
HL = f!(Xc, :L)
@test HL == Y
@test tril(Xc) == tril(Y)
end
end
@test_throws DimensionMismatch hermitianpart(ones(1,2))
for T in (Float64, ComplexF64), uplo in (:U, :L)
A = [randn(T, 2, 2) for _ in 1:2, _ in 1:2]
Aherm = hermitianpart(A, uplo)
@test Aherm == Aherm.data == (A + A')/2
@test Aherm isa Hermitian
@test Aherm.uplo == LinearAlgebra.char_uplo(uplo)
end
end

end # module TestSymmetric

2 comments on commit faf7954

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Executing the daily package evaluation, I will reply here when finished:

@nanosoldier runtests(isdaily = true)

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

Please sign in to comment.