Disambiguate division and pseudo-division?

[martinholters]

...where division means "equivalent to multiplication with inverse" and pseudo-division means "equivalent to multiplication with pseudo-inverse" (with equivalence neglecting any numerical issues). Briefly came up in JuliaLang/julia#23067. Motivating example:

julia> using LinearAlgebra

julia> foo(n) = zeros(n, 3) \ zeros(n, 1)
foo (generic function with 1 method)

julia> foo(1)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

julia> foo(2)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

julia> foo(3)
ERROR: SingularException(1)
Stacktrace:
 [1] ldiv!(::Diagonal{Float64,Array{Float64,1}}, ::Array{Float64,2}) at stdlib/v1.1/LinearAlgebra/src/diagonal.jl:438
 [2] \(::Diagonal{Float64,Array{Float64,1}}, ::Array{Float64,2}) at stdlib/v1.1/LinearAlgebra/src/diagonal.jl:445
 [3] \(::Array{Float64,2}, ::Array{Float64,2}) at stdlib/v1.1/LinearAlgebra/src/generic.jl:862
 [4] foo(::Int64) at ./REPL[2]:1
 [5] top-level scope at none:0

julia> foo(4)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

Here, it would be useful if there was a dedicated operator to ask for the least-squares solution. But if there was one, \ should probably always mean exact solution and fail for rectangular matrices.

So if we do this, it would be quite breaking. Is it justified? Should we explore this?

[yuyichao]

Just pointing out that I complained about this right soon you made this change JuliaLang/julia#23067 (comment) but no one seems to think the ambiguity was an issue for anyone. Also JuliaLang/julia#25431

[Keno]

I wouldn't have any problems with inv(A)*B and pinv(A)*B, where both of those are lazy wrappers. We could even have syntax sugar like A⁻¹B and A⁺B for those.

[andreasnoack]

Here, it would be useful if there was a dedicated operator to ask for the least-squares solution.

This is specified via the the choice of factorization. E.g.

julia> bar(n) = qr(zeros(n, 3), Val(true)) \ zeros(n, 1)
bar (generic function with 1 method)

julia> bar(1)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

julia> bar(2)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

julia> bar(3)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

julia> bar(4)
3×1 Array{Float64,2}:
 0.0
 0.0
 0.0

This should probably be better documented though.

where both of those are lazy wrappers

I'd be reluctant to add more lazy array-like types. It will require a lot of work to support all the combinations. We still have quite some work to do because we are done with Adjoint/Transpose.

[Keno]

I'd be reluctant to add more lazy array-like types. It will require a lot of work to support all the combinations. We still have quite some work to do because we are done with Adjoint/Transpose.

Yeah, though we may have to figure out something more general here.

[goretkin]

I think this ambiguity bit me (for a moment) today. I intended to do a ./ b, but instead wrote a / b, and expected an error for a::Vector, b::Vector.

julia> [1, 1] / [1, 1]
2×2 Array{Float64,2}:
 0.5  0.5
 0.5  0.5

[antoine-levitt]

As @martinholters commented in JuliaLang/julia#44358, it'd be nice to get rid of the DWIM-y thing inherited from Matlab in 2.0. No clear idea of what the syntax for pinvs would like look though. I suggested /+ in another thread, but it clashes with the unary plus. +/ is OK and possibly nice though?