JuliaStats / PDMats.jl

Uniform Interface for positive definite matrices of various structures
Other
104 stars 43 forks source link

Support StaticArrays in X(t)_(inv)A_X(t) with ScalMat #181

Closed devmotion closed 11 months ago

devmotion commented 11 months ago

Currently, ScalMat does not support X_A_Xt etc. with static arrays. For instance, running the tests in this PR on the master branch yields

StaticArrays: Error During Test at REPL[10]:42
  Test threw exception
  Expression: X_A_Xt(A, X) isa SMatrix{10, 10, Float64}
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:43
  Test threw exception
  Expression: X_A_Xt(A, X) ≈ Matrix(X) * Matrix(A) * (Matrix(X))'
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:45
  Test threw exception
  Expression: X_invA_Xt(A, X) isa SMatrix{10, 10, Float64}
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:46
  Test threw exception
  Expression: X_invA_Xt(A, X) ≈ Matrix(X) * (Matrix(A) \ (Matrix(X))')
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:48
  Test threw exception
  Expression: Xt_A_X(A, Y) isa SMatrix{10, 10, Float64}
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:49
  Test threw exception
  Expression: Xt_A_X(A, Y) ≈ (Matrix(Y))' * Matrix(A) * Matrix(Y)
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:51
  Test threw exception
  Expression: Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64}
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
StaticArrays: Error During Test at REPL[10]:52
  Test threw exception
  Expression: Xt_invA_X(A, Y) ≈ (Matrix(Y))' * (Matrix(A) \ Matrix(Y))
  setindex!(::SMatrix{10, 10, Float64, 100}, value, ::Int) is not defined.
   Hint: Use `MArray` or `SizedArray` to create a mutable static array
...
Test Summary: | Pass  Error  Total  Time
StaticArrays  |   54      8     62  4.6s

The PR fixes this issue by changing the implementation to a non-mutating alternative, as done already e.g. for PDiagMat.

The disadvantage (also with the existing code for other psd matrix types) is that this leads to additional allocations when working with e.g. Array. I wonder if generally the non-mutating implementation should be the default (since it works with any matrix type) and additionally possibly internally mutating code for Array should be added (since it avoids allocations and this particular matrix type is so common).

codecov-commenter commented 11 months ago

Codecov Report

All modified lines are covered by tests :white_check_mark:

Files Coverage Δ
src/scalmat.jl 96.82% <100.00%> (ø)

... and 1 file with indirect coverage changes

:loudspeaker: Thoughts on this report? Let us know!.

devmotion commented 11 months ago

The PR is ready for review, @andreasnoack.

I wonder if generally the non-mutating implementation should be the default (since it works with any matrix type) and additionally possibly internally mutating code for Array should be added (since it avoids allocations and this particular matrix type is so common).

That's what I ended up doing for now. I only added specializations for Matrix - originally I thought StridedMatrix might be an alternative but AFAICT it is not guaranteed that a StridedMatrix is mutable.