Skip to content

Solving linear systems with non-trivial RHS #904

@juliohm

Description

@juliohm

As far I understand it, it is always possible to solve a linear system when the entries of the RHS live in a vector space. I have an exotic vector space of compositions and the following happens on Julia v1.6:

julia> using CoDa

julia> b = rand(Composition{3}, 10)
10-element Vector{Composition{3, PARTS} where PARTS}:
 "0.610 : 0.296 : 0.094"
 "0.608 : 0.264 : 0.128"
 "0.510 : 0.016 : 0.473"
 "0.165 : 0.192 : 0.643"
 "0.007 : 0.963 : 0.030"
 "0.698 : 0.255 : 0.048"
 "0.679 : 0.004 : 0.316"
 "0.039 : 0.695 : 0.266"
 "0.763 : 0.097 : 0.140"
 "0.212 : 0.202 : 0.586"

julia> A = rand(10,10)
10×10 Matrix{Float64}:
 0.760859    0.104603   0.823341   0.705046   0.600313    0.916079  0.791862   0.632603    0.778659  0.473457
 0.00990181  0.329258   0.0765397  0.68822    0.536604    0.127633  0.411928   0.659145    0.909937  0.612008
 0.805942    0.498786   0.517988   0.0988711  0.00978311  0.958379  0.322744   0.944924    0.257193  0.970391
 0.212351    0.0373405  0.311988   0.712594   0.30764     0.721235  0.581307   0.00426507  0.681214  0.474252
 0.835093    0.88451    0.832915   0.778178   0.661449    0.139439  0.375395   0.301075    0.289008  0.960027
 0.258823    0.441132   0.268645   0.596426   0.331703    0.675858  0.685603   0.318803    0.812066  0.669362
 0.638684    0.0211262  0.68885    0.0161966  0.309542    0.148414  0.319035   0.866017    0.379569  0.199479
 0.0585245   0.477418   0.886183   0.269309   0.489928    0.750661  0.881117   0.192531    0.234592  0.554885
 0.625465    0.51341    0.157735   0.927946   0.961164    0.390173  0.0871444  0.753436    0.365655  0.100547
 0.102124    0.929157   0.868679   0.234841   0.624787    0.312557  0.0458109  0.440861    0.20748   0.92118

julia> A'A \ A'b # doesn't compile because A'b is a vector of Composition objects (objects in a vector space)
ERROR: MethodError: no method matching one(::Type{Composition})
Closest candidates are:
  one(::Union{Type{T}, T}) where T<:AbstractString at strings/basic.jl:262
  one(::Union{Type{P}, P}) where P<:Dates.Period at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:54
  one(::LinearAlgebra.Tridiagonal{T, V} where V<:AbstractVector{T}) where T at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/special.jl:307
  ...
Stacktrace:
 [1] oneunit(#unused#::Type{Composition})
   @ Base ./number.jl:319
 [2] \(F::LinearAlgebra.LU{Float64, Matrix{Float64}}, B::Vector{Composition})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:101
 [3] \(A::Matrix{Float64}, B::Vector{Composition})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [4] top-level scope
   @ REPL[37]:1

julia> (A'A \ A') * b # manually fix the issue by changing the order of operations
10-element Vector{Composition}:
 "0.000 : 1.000 : 0.000"
 "0.000 : 0.000 : 1.000"
 "0.000 : 0.000 : 1.000"
 "0.000 : 0.000 : 1.000"
 "0.000 : 1.000 : 0.000"
 "0.000 : 1.000 : 0.000"
 "0.000 : 1.000 : 0.000"
 "0.000 : 0.000 : 1.000"
 "0.992 : 0.008 : 0.000"
 "0.000 : 1.000 : 0.000"

Would it make sense to generalize the backslash operator to accept non-trivial RHS? The workaround above is more expensive.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions