Improve StaticArray + BandedMatrices support (update of #131) (#152)

* Don't promote to `Matrix`.

* Add StaticArrays tests, fix A / B::PDiagMat

* Improve `StaticArray` support and add tests

* Bump version

* Fix merge errors

* Bump version

* Add tests with BandedMatrices

* Fix test error with old Julia versions

Co-authored-by: chriselrod <[email protected]>

Author: devmotion

Project.toml

@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.6"
+version = "0.11.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -11,7 +11,9 @@ SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 julia = "1"
 
 [extras]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["BandedMatrices", "StaticArrays", "Test"]

src/pdiagmat.jl

@@ -139,25 +139,25 @@ end
 
 ### tri products
 
-function X_A_Xt(a::PDiagMat, x::StridedMatrix)
+function X_A_Xt(a::PDiagMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 2)
     z = x .* sqrt.(permutedims(a.diag))
     z * transpose(z)
 end
 
-function Xt_A_X(a::PDiagMat, x::StridedMatrix)
+function Xt_A_X(a::PDiagMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 1)
     z = x .* sqrt.(a.diag)
     transpose(z) * z
 end
 
-function X_invA_Xt(a::PDiagMat, x::StridedMatrix)
+function X_invA_Xt(a::PDiagMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 2)
     z = x ./ sqrt.(permutedims(a.diag))
     z * transpose(z)
 end
 
-function Xt_invA_X(a::PDiagMat, x::StridedMatrix)
+function Xt_invA_X(a::PDiagMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 1)
     z = x ./ sqrt.(a.diag)
     transpose(z) * z

src/pdmat.jl

@@ -16,8 +16,7 @@ function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S}
     PDMat{T,S}(d, convert(S, mat), chol)
 end
 
-PDMat(mat::Matrix) = PDMat(mat, cholesky(mat))
-PDMat(mat::Symmetric) = PDMat(Matrix(mat))
+PDMat(mat::AbstractMatrix) = PDMat(mat, cholesky(mat))
 PDMat(fac::Cholesky) = PDMat(Matrix(fac), fac)
 
 ### Conversion
@@ -94,25 +93,25 @@ invquad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.ma
 
 ### tri products
 
-function X_A_Xt(a::PDMat, x::StridedMatrix)
+function X_A_Xt(a::PDMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 2)
     z = x * chol_lower(a.chol)
     return z * transpose(z)
 end
 
-function Xt_A_X(a::PDMat, x::StridedMatrix)
+function Xt_A_X(a::PDMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 1)
     z = chol_upper(a.chol) * x
     return transpose(z) * z
 end
 
-function X_invA_Xt(a::PDMat, x::StridedMatrix)
+function X_invA_Xt(a::PDMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 2)
     z = x / chol_upper(a.chol)
     return z * transpose(z)
 end
 
-function Xt_invA_X(a::PDMat, x::StridedMatrix)
+function Xt_invA_X(a::PDMat, x::AbstractMatrix)
     @check_argdims dim(a) == size(x, 1)
     z = chol_lower(a.chol) \ x
     return transpose(z) * z

test/runtests.jl

@@ -1,5 +1,5 @@
 include("testutils.jl")
-tests = ["pdmtypes", "addition", "generics", "kron", "chol"]
+tests = ["pdmtypes", "addition", "generics", "kron", "chol", "specialarrays"]
 println("Running tests ...")
 
 for t in tests

test/specialarrays.jl

@@ -0,0 +1,72 @@
+using BandedMatrices
+using StaticArrays
+
+@testset "Special matrix types" begin
+    @testset "StaticArrays" begin
+        # Full matrix
+        S = (x -> x * x')(@SMatrix(randn(4, 7)))
+        PDS = PDMat(S)
+        @test PDS isa PDMat{Float64, <:SMatrix{4, 4, Float64}}
+        @test isbits(PDS)
+
+        # Diagonal matrix
+        D = PDiagMat(@SVector(rand(4)))
+        @test D isa PDiagMat{Float64, <:SVector{4, Float64}}
+
+        x = @SVector rand(4)
+        X = @SMatrix rand(10, 4)
+        Y = @SMatrix rand(4, 10)
+
+        for A in (PDS, D)
+            @test A * x isa SVector{4, Float64}
+            @test A * x ≈ Matrix(A) * Vector(x)
+
+            @test A * Y isa SMatrix{4, 10, Float64}
+            @test A * Y ≈ Matrix(A) * Matrix(Y)
+
+            @test X / A isa SMatrix{10, 4, Float64}
+            @test X / A ≈ Matrix(X) / Matrix(A)
+
+            @test A \ x isa SVector{4, Float64}
+            @test A \ x ≈ Matrix(A) \ Vector(x)
+
+            @test A \ Y isa SMatrix{4, 10, Float64}
+            @test A \ Y ≈ Matrix(A) \ Matrix(Y)
+
+            @test X_A_Xt(A, X) isa SMatrix{10, 10, Float64}
+            @test X_A_Xt(A, X) ≈ Matrix(X) * Matrix(A) *  Matrix(X)'
+
+            @test X_invA_Xt(A, X) isa SMatrix{10, 10, Float64}
+            @test X_invA_Xt(A, X) ≈ Matrix(X) * (Matrix(A) \ Matrix(X)')
+
+            @test Xt_A_X(A, Y) isa SMatrix{10, 10, Float64}
+            @test Xt_A_X(A, Y) ≈ Matrix(Y)' * Matrix(A) * Matrix(Y)
+
+            @test Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64}
+            @test Xt_invA_X(A, Y) ≈ Matrix(Y)' * (Matrix(A) \ Matrix(Y))
+        end
+    end
+
+    @testset "BandedMatrices" begin
+        # Full matrix
+        A = Symmetric(BandedMatrix(Eye(5), (1, 1)))
+        P = PDMat(A)
+        @test P isa PDMat{Float64, <:BandedMatrix{Float64}}
+
+        x = rand(5)
+        X = rand(2, 5)
+        Y = rand(5, 2)
+        @test P * x ≈ A * x
+        @test P * Y ≈ A * Y
+        # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594
+        if VERSION >= v"1.3.0-DEV.562"
+            @test X / P ≈ X / A
+        end
+        @test P \ x ≈ A \ x
+        @test P \ Y ≈ A \ Y
+        @test X_A_Xt(P, X) ≈ X * A * X'
+        @test X_invA_Xt(P, X) ≈ X * (A \ X')
+        @test Xt_A_X(P, Y) ≈ Y' * A * Y
+        @test Xt_invA_X(P, Y) ≈ Y' * (A \ Y)
+    end
+end