Underdetermined Linear System General Solution

Specific Domains Numerics
1

Hey all, this may be more of a question about linear algebra than it is about Julia. Given an underdetermined linear system A x = b, where A has dimensions (m, n) we expect there to be a vector space of dimension (m - n) describing all x that solve our system (right?).

How do I get the basis vectors for that subspace? I know A \ b will give me a particular x, and I could then try orthogonal perturbations to see if they are still valid solutions, but surely there is a less ugly way.

Thanks!

2

So you have Ax=b where A\in \mathbb{R}^{m\times n} with m<n? Then let x_1 be one particular solution, say, found by x1 = A\b. Then, any other value x_1 + x_0 is also a solution if x_0 solves Ax_0=0. All possible vectors x_0 solving Ax_0=0 are linear combinations of basis vectors in the nullspace of A, i.e., x_0 \in {\cal N}(A).
Example:

julia> using LinearAlgebra

julia> A = rand(-9:9,2,3)
2×3 Array{Int64,2}:
 -7  -9  -7
 -8  -8  -1

julia> rank(A)
2

julia> b = [1,2]
2-element Array{Int64,1}:
 1
 2

julia> x1 = A\b
3-element Array{Float64,1}:
 -0.1637895602137279
 -0.10583641594739004
  0.1570078092889437

julia> N = nullspace(A)
3×1 Array{Float64,2}:
 -0.673770210645276
  0.7024412834386926
 -0.22936858234732868

julia> x0 = 0.5*N[:]
3-element Array{Float64,1}:
 -0.336885105322638
  0.3512206417193463
 -0.11468429117366434

julia> A*(x1+x0)
2-element Array{Float64,1}:
 0.9999999999999996
 1.9999999999999971

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

So… by doing N = nullspace(A), you find a matrix N containing basis vectors for the nullspace (kernel) of A. Then, with x1 = A\b, x1 is one solution to Ax=b, but x1+N*z is also a solution to Ax+b for any vector z.

6 Likes
3

Thank you @BLI!

4

I could add that A\b gives the solution x with minimal norm (uses Singular Value Decomposition, I think). If you find x by QR decomposition, you can choose a solution with the most zeros in x if that is more preferable.

1 Like
5

The minimum norm solutions is computed with the pivoted QR followed by a complete orthogonal factorization of the trapezoid R factor, see Complete Orthogonal Factorization.

3 Likes
6

Would you mind sharing how you’d do that?

7

Well, I guess QR factorization of A\in \mathbb{R}^{m\times n} with m<n leads to A=QR=b \Rightarrow R=Q'b = \tilde{b} since Q^{-1} = Q' (Q' is the transpose of Q).

Because A has fewer rows than columns and R is upper triangular, R has structure R = (R_1, M) where R_1 is upper triangular and square, while M may be anything.

With R_1 upper triangular, it has rank m, and we can split x into x_1 of length m and x_2 of length n-m.

Then Ax=b \Rightarrow R_1 x_1 + M x_2 = \tilde{b}, or x_1 = R_1^{-1}(\tilde{b}-Mx_2)

Here, you can choose x_2 to be anything you want, e.g., x_2=0, and find a unique x_1 which satisfies Ax=b. With this procedure, x will be guaranteed to hold at least n-m zeros.

A simple example:

julia> using LinearAlgebra

julia> A = rand(-9:9,2,4)
2×4 Array{Int64,2}:
 -1   2   5  -2
 -6  -4  -1   5

julia> b=rand(-9:9,2)
2-element Array{Int64,1}:
 -5
 -8

julia> rank(A)
2

julia> F = qr(A);
julia> F.R
2×4 Array{Float64,2}:
 6.08276   3.61678   0.164399  -4.60317
 0.0      -2.63038  -5.09637    2.79478

julia> R1 = F.R[:,1:2]
2×2 Array{Float64,2}:
 6.08276   3.61678
 0.0      -2.63038

julia> M = F.R[:,3:end]
2×2 Array{Float64,2}:
  0.164399  -4.60317
 -5.09637    2.79478

julia> x1 = R1\F.Q'*b
2-element Array{Float64,1}:
  2.2500000000000027
 -1.3750000000000024

julia> x = vcat(x1,zeros(2))
4-element Array{Float64,1}:
  2.2500000000000027
 -1.3750000000000024
  0.0
  0.0

julia> A*x
2-element Array{Float64,1}:
 -5.000000000000007
 -8.000000000000007

julia> b
2-element Array{Int64,1}:
 -5
 -8

Here, I haven’t had use of matrix M — since I chose x_2=0. However, if you for some reason is a fan of the Hitchhikers Guide to the Galaxy and like the number 42, you may do as follows:

julia> x2 = fill(42,2)
2-element Array{Int64,1}:
 42
 42

julia> x1 = R1\(F.Q'*b-M*x2)
2-element Array{Float64,1}:
  54.75000000000007
 -38.125000000000064

julia> x = vcat(x1,x2)
4-element Array{Float64,1}:
  54.75000000000007
 -38.125000000000064
  42.0
  42.0

julia> A*x
2-element Array{Float64,1}:
 -5.000000000000199
 -8.000000000000199

julia> b
2-element Array{Int64,1}:
 -5
 -8

NOTE: it is possible that \mathrm{rank}(A) < \min(m,n). In that case, you need to split up R further.

2 Likes
8

Note however that that only happens for dense arrays, which is a bit of a catch, at least among students:

julia> using  LinearAlgebra, SparseArrays

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

julia> norm(A \ b)
0.5187887610670437

julia> norm(sparse(A) \ b)
2.433310933552997

For sparse matrices, the QR factorization is performed by SuiteSparseQR (an external C library), which returns a basic solution instead of a minimum-norm solution. A basic solution is one with many zeros (and often does not have minimum norm).

3 Likes