StaticArrays solve symmetric linear system - seems typeinstable

General Usage Performance
1

Hi there, I am wondering why the following does not return a SVector, but allocates a normal Vector. Also warntype shows problems

using StaticArrays
using LinearAlgebra
using Statistics
function testme(xs=rand(20))
    ys = xs.*4
    A = Symmetric(SA_F64[
        mean(x*x for x in xs)  mean(xs)
        mean(xs)               1.0
    ])
    b = SA_F64[mean(x*y for (x, y) in zip(xs, ys)), mean(ys)]
    A \ b
end
typeof(testme())  # Vector, not SVector
@code_warntype testme(rand(20))

here the output of warntype, where the last ::Union{Vector{Float64}, SVector{2, Float64}} is highlighted in yellow

1 ─ %1  = Base.broadcasted(Main.:*, xs, 4)::Core.PartialStruct(Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(*), Tuple{Vector{Float64}, Int64}}, Any[Core.Const(Base.Broadcast.DefaultArrayStyle{1}()), Core.Const(*), Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(4)]), Nothing])
│         (ys = Base.materialize(%1))
│   %3  = Main.SA_F64::Core.Const(SA_F64)
│   %4  = Core.tuple(2, 2)::Core.Const((2, 2))
│         (#31 = %new(Main.:(var"#31#33")))
%6  = #31::Core.Const(var"#31#33"())=======>            ]  30/43
│   %7  = Base.Generator(%6, xs)::Base.Generator{Vector{Float64}, var"#31#33"}
│   %8  = Main.mean(%7)::Float64
│   %9  = Main.mean(xs)::Float64BaseExt
│   %10 = Main.mean(xs)::Float64t
│   %11 = Base.typed_hvcat(%3, %4, %8, %9, %10, 1.0)::Core.PartialStruct(SMatrix{2, 2, Float64, 4}, Any[Core.PartialStruct(NTuple{4, Float64}, Any[Float64, Float64, Float64, Core.Const(1.0)])])rsMakieExt
│         (A = Main.Symmetric(%11))
│   %13 = Main.SA_F64::Core.Const(SA_F64)
│         (#32 = %new(Main.:(var"#32#34")))
│   %15 = #32::Core.Const(var"#32#34"())
│   %16 = Main.zip(xs, ys)::Base.Iterators.Zip{Tuple{Vector{Float64}, Vector{Float64}}}
│   %17 = Base.Generator(%15, %16)::Base.Generator{Base.Iterators.Zip{Tuple{Vector{Float64}, Vector{Float64}}}, var"#32#34"}
│   %18 = Main.mean(%17)::Float64
│   %19 = Main.mean(ys)::Float64
│         (b = Base.getindex(%13, %18, %19))
│   %21 = (A::Core.PartialStruct(Symmetric{Float64, SMatrix{2, 2, Float64, 4}}, Any[Core.PartialStruct(SMatrix{2, 2, Float64, 4}, Any[Core.PartialStruct(NTuple{4, Float64}, Any[Float64, Float64, Float64, Core.Const(1.0)])]), Core.Const('U')]) \ b)::Union{Vector{Float64}, SVector{2, Float64}}
└──       return %21

Without Symmetric, everything is fine and gives a SVector correctly.
Any help is highly appreciated

2

The problem is the poly-algorithm of \ for Symmetric ends up using a Bunch–Kaufmann factorization, which isn’t implemented in StaticArrays and falls back to the generic LinearAlgebra routine that allocates a Vector for the result. Probably you want

cholesky(A) \ b

since your A in this case is positive-definite (because your A = X^T X / n is from a least-square problem). Cholesky factorization is supported by StaticArrays.

3 Likes
3

thank you a lot for this compact deep dive

1 Like