Correct isposdef for CompositeMap (#151)

Author: JeffFessler

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.3.1"
+version = "3.4.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/composition.jl

@@ -21,8 +21,7 @@ Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessar
 
 # the following rules are sufficient but not necessary
 for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
-                    (:ishermitian, :_ishermitian, :adjoint),
-                    (:isposdef, :_isposdef, :adjoint))
+                    (:ishermitian, :_ishermitian, :adjoint))
     @eval begin
         LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
         $_f(maps::Tuple{}) = true
@@ -42,6 +41,15 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
     end
 end
 
+# A * B * A and A * B * A' are positive definite if (sufficient condition) A & B are positive definite
+LinearAlgebra.isposdef(A::CompositeMap) = _isposdef(A.maps)
+_isposdef(maps::Tuple{}) = true # empty product is equivalent to "I" which is pos. def.
+_isposdef(maps::Tuple{<:LinearMap}) = isposdef(maps[1])
+function _isposdef(maps::Tuple{Vararg{LinearMap}})
+    (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) && 
+    isposdef(maps[1]) && _isposdef(Base.front(Base.tail(maps)))
+end
+
 # scalar multiplication and division (non-commutative case)
 function Base.:(*)(α::Number, A::LinearMap)
     T = promote_type(typeof(α), eltype(A))

test/composition.jl

@@ -28,7 +28,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     @test @inferred F4 * v == @inferred F * (F * (F * (F * v)))
     @test @inferred Matrix(M * transpose(M)) ≈ A * transpose(A)
     @test @inferred !isposdef(M * transpose(M))
-    @test @inferred isposdef(M * M')
+    @test @inferred isposdef(LinearMap(M * M', isposdef=true))
     @test @inferred issymmetric(N * N')
     @test @inferred ishermitian(N * N')
     @test @inferred !issymmetric(M' * M)
@@ -40,10 +40,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     @test @inferred !issymmetric(FC'FC)
     @test @inferred ishermitian(FC'FC)
     @test @inferred ishermitian(FC'*H*FC)
-    @test @inferred isposdef(transpose(F) * F * 3)
-    @test @inferred isposdef(transpose(F) * 3 * F)
+    @test @inferred issymmetric(transpose(F) * F * 3)
+    @test @inferred issymmetric(transpose(F) * 3 * F)
     @test @inferred !isposdef(-5*transpose(F) * F)
-    @test @inferred isposdef((M * F)' * M * 4 * F)
+    @test @inferred ishermitian((M * F)' * M * 4 * F)
     @test @inferred transpose(M * F) == @inferred transpose(F) * transpose(M)
     @test @inferred (4*((-3*M)*2)) == @inferred -12M*2
     @test @inferred (4*((3*(-M))*2)*(-5)) == @inferred -12M*(-10)
@@ -132,4 +132,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     X = rand(size(Ls, 2), 10)
     Y = similar(X, (size(Ls, 1), size(X, 2)))
     @test mul!(Y, Ls, X) ≈ L4.lmap * L3.lmap * L2.lmap * L1.lmap * X
+
+    # test isposdef on a case where sufficient conditions work
+    B = LinearMap([1 0; 0 1], isposdef=true) # isposdef!
+    C = B' * B * B * B * B # no B' at end on purpose
+    @test @inferred isposdef(C)
+    @test @inferred isposdef(B * B) # even case for empty tuple test
 end

test/scaledmap.jl

@@ -77,9 +77,10 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     @test mul!(x2, C', y1) == x1
 
     # check scale*conj(scale)
-    A = LinearMap{Float32}(rand(N, 2)) # rank=2 w.p.1
+    A = LinearMap{Float32}(rand(1, 2)) # rank=1 w.p.1
     B = @inferred β * A
     C = @inferred B' * B
+    @test !isposdef(C)
     @test @inferred isposdef(C.λ)
 
     N = 2^8