How can I use \= with factorizations?

General Usage
1

As per title.

For example, this works:

using LinearAlgebra
A = rand(3, 3)
b = rand(3)
@show b
b \= A
@show b

b = [0.09019616214492654, 0.2564457656075061, 0.7120782534101571]
b = [0.5468125936523385 1.4817879951917017 0.5732261599468792]

but this doesn’t:

using LinearAlgebra
A = factorize(rand(3, 3))
b = rand(3)
@show b
b \= A
@show b

b = [0.54431634700671, 0.6989662716060527, 0.244305881937394]
StackOverflowError:

Stacktrace:
     [1] \(x::Vector{Float64}, y::LU{Float64, Matrix{Float64}})
       @ Base ./operators.jl:598
     [2] /(A::Adjoint{Float64, LU{Float64, Matrix{Float64}}}, B::Adjoint{Float64, Vector{Float64}})
       @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1144
--- the last 2 lines are repeated 25974 more times ---
 [51951] \(x::Vector{Float64}, y::LU{Float64, Matrix{Float64}})
       @ Base ./operators.jl:598
 [51952] top-level scope
       @ In[69]:5
 [51953] eval
       @ ./boot.jl:360 [inlined]

What’s the easiest way to make the second one work too? Should I add a new method to ldiv!? Would that be enough?

2

It gives an answer, but I’m not sure it does what you think it does.

b \= A is equivalent to b = b \ A. This means two things:

  1. It is essentially equivalent to pinv(b) * A — that is, you are dividing A by b, not b by A. Did you want b' / A instead, i.e. b^T A^{-1}?

  2. It is not in-place. It first computes b \ A, allocating a new row vector for the result, and then assigns the result to b (discarding the original column vector).

3

Are you sure about (1.)? pinv(b) * A gives me a DimensionMismatch("A has dimensions (3,1) but B has dimensions (3,3)"). I assumed that was equivalent to b = A \ b, which dimensionally makes sense.

4

The docstring says you’re wrong, so if it does something else that’s probably a bug.

5

Quite sure:

julia> A = rand(3, 3); b = rand(3);

julia> pinv(b)*A
1×3 transpose(::Vector{Float64}) with eltype Float64:
 0.998627  1.23721  1.55637

julia> b \ A
1×3 transpose(::Vector{Float64}) with eltype Float64:
 0.998627  1.23721  1.55637

I think you were running this after you ran b \= A, which replaced b (a column vector) with the row vector b \ A. And using a row vector, e.g. pinv(b')*A, indeed gives a dimension-mismatch error, because this doesn’t work with the definition of the pseudo-inverse.

1 Like
6

This is never true in Julia. x OP= y is always equivalent to x = x OP y.

7

Ah! I see. Then the question comes down to: Is there a way to add a new definition for \=?

8

No. x \= y is just syntactic sugar for x = x \ y. There is no way to redefine \= by itself.

1 Like