How can I use \= with factorizations?

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
@ 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?

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).

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.

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

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

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

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

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

1 Like