Add matrix symmetrizing functions

This adds functions `symmetrize` and `symmetrize!`, which produce
symmetric matrices based on the input `X` via efficient computation of
`(X + X') / 2`.

Author: ararslan

NEWS.md

@@ -29,6 +29,7 @@ Standard library changes
 #### LinearAlgebra
 
 * `diagm` and `spdiagm` now accept optional `m,n` initial arguments to specify a size ([#31654]).
+* New functions `symmetrize` and `symmetrize!` for constructing symmetric matrices ([#TODO]).
 
 #### SparseArrays
 

stdlib/LinearAlgebra/docs/src/index.md

@@ -430,6 +430,8 @@ Base.copy(::Union{Transpose,Adjoint})
 LinearAlgebra.stride1
 LinearAlgebra.checksquare
 LinearAlgebra.peakflops
+LinearAlgebra.symmetrize
+LinearAlgebra.symmetrize!
 ```
 
 ## Low-level matrix operations

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -130,6 +130,8 @@ export
     svdvals!,
     svdvals,
     sylvester,
+    symmetrize,
+    symmetrize!,
     tr,
     transpose,
     transpose!,

stdlib/LinearAlgebra/src/symmetric.jl

@@ -890,3 +890,38 @@ end
 *(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B
 *(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractTriangular}) = A.parent * B
 *(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Transpose{<:Any,<:AbstractTriangular}) = A.parent * B
+
+"""
+    symmetrize!(X::AbstractMatrix{<:AbstractFloat})
+
+Make the matrix `X` symmetric in-place using the formula `(X + X') / 2`.
+
+See also: [`symmetrize`](@ref)
+
+!!! compat "Julia 1.3"
+    This function requires Julia 1.3 or later.
+"""
+function symmetrize!(X::AbstractMatrix{<:AbstractFloat})
+    require_one_based_indexing(X)
+    n = checksquare(X)
+    @inbounds for j = 2:n, i = 1:j-1
+        X[i,j] = (X[i,j] + X[j,i]) / 2
+        X[j,i] = X[i,j]
+    end
+    X
+end
+
+"""
+    symmetrize(X::AbstractMatrix)
+
+Construct a symmetric matrix based on `X` using the formula `(X + X') / 2`.
+This function is a no-op for [`Symmetric`](@ref)-wrapped matrices.
+
+See also: [`symmetrize!`](@ref)
+
+!!! compat "Julia 1.3"
+    This function requires Julia 1.3 or later.
+"""
+symmetrize(X::AbstractMatrix{<:AbstractFloat}) = symmetrize!(copy(X))
+symmetrize(X::AbstractMatrix{<:Real}) = symmetrize!(float(X))
+symmetrize(X::Symmetric) = X

stdlib/LinearAlgebra/test/symmetric.jl

@@ -603,4 +603,26 @@ end
     @test LinearAlgebra.hermitian_type(Int) == Int
 end
 
+@testset "symmetrize" begin
+    @testset "in-place" begin
+        X = rand(Float32, 4, 4)
+        Y = copy(X)
+        symmetrize!(X)
+        @test X ≈ (Y + Y') / 2
+        @test issymmetric(X)
+    end
+    @testset "out-of-place" begin
+        X = randn(Float32, 4, 4)
+        Y = symmetrize(X)
+        @test Y ≈ (X + X') / 2
+        @test issymmetric(Y)
+        X = Int32[1 2 3; 4 5 6; 7 8 9]
+        @test symmetrize(X) ≈ Float32[1 3 5; 3 5 7; 5 7 9]
+    end
+    @testset "Symmetric" begin
+        S = Symmetric([1 2; 2 1])
+        @test symmetrize(S) === S
+    end
+end
+
 end # module TestSymmetric