Matrix division `\` vs compute by hand, why such big difference?

General Usage
1

In the following example, it seems that the Matrix division \ does not provide an intended computation.

Let x= A'\b, then if norm(A'*x-b, Inf) is zero such that A'x-b == 0, then b is a linear combination of rows from A (b lies in the row space of A).

I am using Matrix division \ to test if b lies in the row space of A or not, and find that the \ gives a result different from the result by the reduced row echelon form. Hence, I compute it by hand, the norm(A'*x-b, Inf) is 4.755666482359269e-12, but the result from Matrix division is 0.00010237, far from zero (tolerance)

Download matH.jls first.

using LinearAlgebra, Serialization


function loadData(fn)
    #load data
    local t
    open(fn, "r") do io
        t = deserialize(io)
    end
    return t
end



function main()
    H, v = loadData("./matH.jls")

    G = qr(H', ColumnNorm())
    z = G.P*inv(G.R)*Matrix(G.Q)'*v
    display(norm(H'*z-v, Inf))  #4.755666482359269e-12

    z1 = G\v
    display(norm(H'*z1-v, Inf)) #0.00010237343843657264

end

main()

nothing

From the \(A,B) , the docstring of Julia: " For rectangular A the result is the minimum-norm least squares solution computed by a pivoted QR factorization of A and a rank estimate of A based on the R factor."

Is there any setting I can choose when using Matrix division \ in Julia?

Edit: more discussions on `x=A\b` by `qr(A, ColumnNorm()) \ B`, Matrix division does not yield `Ax=b` · Issue #49625 · JuliaLang/julia · GitHub

1 Like