Speed up WrappedMap constructor and LM-promotion (#154)

Author: dkarrasch

docs/src/history.md

@@ -1,5 +1,16 @@
 # Version history
 
+## What's new in v3.4
+
+* In `WrappedMap` constructors, as implicitly called in addition and mutliplication
+  of `LinearMap`s and `AbstractMatrix` objects, (conjugate) symmetry and positive
+  definiteness are only determined for matrix types for which these checks are expected
+  to be very cheap or even known at compile time based on the concrete type. The default
+  for `LinearMap` subtypes is to call, for instance, `issymmetric`, because symmetry
+  properties are either stored or easily obtained from constituting maps. For custom matrix
+  types, define corresponding methods `LinearMaps._issymmetric`, `LinearMaps._ishermitian`
+  and `LinearMaps._isposdef` to hook into the property checking mechanism.
+
 ## What's new in v3.3
 
 * `AbstractVector`s can now be wrapped by a `LinearMap` just like `AbstractMatrix``

src/LinearMaps.jl

@@ -48,7 +48,7 @@ LinearAlgebra.isposdef(::LinearMap) = false # default assumptions
 
 Base.ndims(::LinearMap) = 2
 Base.size(A::LinearMap, n) =
-    (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
+    (n == 1 || n == 2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
 # check dimension consistency for multiplication A*B
@@ -267,24 +267,27 @@ include("show.jl") # show methods for LinearMap objects
 
 Construct a linear map object, either from an existing `LinearMap` or `AbstractVecOrMat` `A`,
 with the purpose of redefining its properties via the keyword arguments `kwargs`;
-a `UniformScaling` object `J` with specified (square) dimension `M`; from a `Number`
-object to lazily represent filled matrices; or
+a `UniformScaling` object `J` with specified (square) dimension `M`; or
 from a function or callable object `f`. In the latter case, one also needs to specify
 the size of the equivalent matrix representation `(M, N)`, i.e., for functions `f` acting
 on length `N` vectors and producing length `M` vectors (with default value `N=M`).
 Preferably, also the `eltype` `T` of the corresponding matrix representation needs to be
-specified, i.e. whether the action of `f` on a vector will be similar to, e.g., multiplying
+specified, i.e., whether the action of `f` on a vector will be similar to, e.g., multiplying
 by numbers of type `T`. If not specified, the devault value `T=Float64` will be assumed.
 Optionally, a corresponding function `fc` can be specified that implements the adjoint
 (=transpose in the real case) of `f`.
 
 The keyword arguments and their default values for the function-based constructor are:
-*   `issymmetric::Bool = false` : whether `A` or `f` acts as a symmetric matrix
-*   `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` acts as a Hermitian
+*   `issymmetric::Bool = false` : whether `A` or `f` act as a symmetric matrix
+*   `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` act as a Hermitian
     matrix
-*   `isposdef::Bool = false` : whether `A` or `f` acts as a positive definite matrix.
+*   `isposdef::Bool = false` : whether `A` or `f` act as a positive definite matrix.
 For existing linear maps or matrices `A`, the default values will be taken by calling
-`issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.
+internal functions `_issymmetric`, `_ishermitian` and `_isposdef` on the existing object `A`.
+These in turn dispatch to (overloads of) `LinearAlgebra`'s `issymmetric`, `ishermitian`,
+and `isposdef` methods whenever these checks are expected to be computationally cheap or even
+known at compile time as for certain structured matrices, but return `false` for generic
+`AbstractMatrix` types.
 
 For the function-based constructor, there is one more keyword argument:
 *   `ismutating::Bool` : flags whether the function acts as a mutating matrix multiplication

src/wrappedmap.jl

@@ -5,12 +5,11 @@ struct WrappedMap{T, A<:MapOrVecOrMat} <: LinearMap{T}
     _isposdef::Bool
 end
 function WrappedMap{T}(lmap::MapOrMatrix;
-                        issymmetric::Bool = issymmetric(lmap),
-                        ishermitian::Bool = ishermitian(lmap),
-                        isposdef::Bool = isposdef(lmap)) where {T}
+                        issymmetric::Bool = _issymmetric(lmap),
+                        ishermitian::Bool = _ishermitian(lmap),
+                        isposdef::Bool = _isposdef(lmap)) where {T}
     WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
 end
-WrappedMap(lmap::MapOrMatrix{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...)
 function WrappedMap{T}(lmap::AbstractVector;
                         issym::Bool = false,
                         isherm::Bool = false,
@@ -20,7 +19,23 @@ function WrappedMap{T}(lmap::AbstractVector;
                                 length(lmap) == 1 && ishermitian(first(lmap)),
                                 length(lmap) == 1 && isposdef(first(lmap)))
 end
-WrappedMap(lmap::AbstractVector{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...)
+WrappedMap(lmap::MapOrVecOrMat{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...)
+
+# cheap property checks (usually by type)
+_issymmetric(A::AbstractMatrix) = false
+_issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
+_issymmetric(A::LinearMap) = issymmetric(A)
+_issymmetric(A::LinearAlgebra.RealHermSymComplexSym) = issymmetric(A)
+_issymmetric(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = issymmetric(A)
+
+_ishermitian(A::AbstractMatrix) = false
+_ishermitian(A::AbstractSparseMatrix) = ishermitian(A)
+_ishermitian(A::LinearMap) = ishermitian(A)
+_ishermitian(A::LinearAlgebra.RealHermSymComplexHerm) = ishermitian(A)
+_ishermitian(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = ishermitian(A)
+
+_isposdef(A::AbstractMatrix) = false
+_isposdef(A::LinearMap) = isposdef(A)
 
 const VecOrMatMap{T} = WrappedMap{T,<:AbstractVecOrMat}
 

test/wrappedmap.jl

@@ -3,13 +3,13 @@ using Test, LinearMaps, LinearAlgebra
 @testset "wrapped maps" begin
     A = rand(10, 20)
     B = rand(ComplexF64, 10, 20)
-    SA = A'A + I
-    SB = B'B + I
+    SA = Hermitian(A'A + I)
+    SB = Hermitian(B'B + I)
     L = @inferred LinearMap{Float64}(A)
     @test summary(L) == "10×20 LinearMaps.WrappedMap{Float64}"
     @test occursin("10×20 LinearMaps.WrappedMap{Float64}", sprint((t, s) -> show(t, "text/plain", s), L))
-    MA = @inferred LinearMap(SA)
-    MB = @inferred LinearMap(SB)
+    MA = @inferred LinearMap(SA, isposdef=true)
+    MB = @inferred LinearMap(SB, isposdef=true)
     @test eltype(Matrix{Complex{Float32}}(LinearMap(A))) <: Complex
     @test size(L) == size(A)
     @test @inferred !issymmetric(L)
@@ -36,7 +36,7 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred LinearMap(M')' * v == A * v
     @test @inferred(transpose(transpose(M))) === M
     @test @inferred(adjoint(adjoint(M))) === M
-    Mherm = @inferred LinearMap(A'A)
+    Mherm = @inferred LinearMap(Hermitian(A'A), isposdef=true)
     @test @inferred ishermitian(Mherm)
     @test @inferred !issymmetric(Mherm)
     @test @inferred !issymmetric(transpose(Mherm))
@@ -96,4 +96,21 @@ using Test, LinearMaps, LinearAlgebra
     @test issymmetric(U)
     @test !ishermitian(U)
     @test !isposdef(U)
+    # symmetry
+    O4 = ones(4)
+    O3 = ones(3)
+    Z3 = zeros(3)
+    for M in (Bidiagonal(O4, O3, :U),
+              Bidiagonal(O4, Z3, :L),
+              Diagonal(O4),
+              Diagonal(im*O4),
+              SymTridiagonal(O4, O3),
+              Tridiagonal(O3, O4, O3),
+              Tridiagonal(O3, O4, Z3),
+              sparse(Tridiagonal(O3, O4, O3)),
+              sparse(Diagonal(im*O4)),
+            )
+        @test issymmetric(LinearMap(M)) == issymmetric(M)
+        @test ishermitian(LinearMap(M)) == ishermitian(M)
+    end
 end