Unexpected solution for fewer equations than unknowns


The following system has many solutions, per simple Gaussian reduction:

x - 3y + z = 1
x + y + 2z = 14

Yet Julia actually gives me a unique solution:

julia> A = [1 -3 1; 1 1 2]; u = [1, 14];

julia> A
2×3 Array{Int64,2}:
 1  -3  1
 1   1  2

julia> u
2-element Array{Int64,1}:

julia> A \ u
3-element Array{Float64,1}:

numpy errors out, as expected.

Am I entering the vectors incorrectly? What’s going on?

You’re getting the minimum-norm solution.

1 Like

Oh. Thank you for pointing that out. This is my first time back in linear algebra in 20 years, I guess I haven’t gotten that far yet.

The left-division doc didn’t mention that at all. I think I’d prefer an error.

Thanks again!

I’m not sure if/where that is explained in the documentation. It’s a side effect of calling the underlying LAPACK functions. When A is rectangular, \ performs a QR factorization. If A has more columns than rows, you receive the minimum-norm solution (in the Euclidean norm). If A has more rows than columns and has full column rank, you receive the least-squares solution. If A has more rows than columns but is column-rank deficient, you receive the minimum-norm least-squares solution. That is often the desired behavior in numerical computations.

However, that behavior goes away when you work with sparse matrices, which is an unfortunate inconsistency:

julia> using SparseArrays

julia> A = [1. -3 1; 1 1 2]; u = [1, 14];  # A must have real elements

julia> B = sparse(A);

julia> B \ u
3-element Array{Float64,1}:

In this case, you receive a basic solution, i.e., one with as many zeros as the column-rank deficiency of A. It’s not the minimum-norm solution.

1 Like

@dpo thanks for the pointer to LAPACK. The main thing I am taking away from this is that what I am seeing is “expected behavior” and the details are best left as an exercise. I am ok with that.