# Underdetermined Linear System General Solution

**URL:** <https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865>\
**Category:** Numerics\
**Created:** [April 19, 2020, 6:20pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865 "2020-04-19T18:20:30Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Alex\_Ellison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex_ellison/32/5749_2.png) [@Alex\_Ellison](https://discourse.julialang.org/u/Alex_Ellison)\
**Post date:** [April 19, 2020, 6:20pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/1 "2020-04-19T18:20:30Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 19, 2020, 6:47pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/2 "2020-04-19T18:47:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Alex\_Ellison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex_ellison/32/5749_2.png) [@Alex\_Ellison](https://discourse.julialang.org/u/Alex_Ellison)\
**Post date:** [April 19, 2020, 7:08pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/3 "2020-04-19T19:08:16Z")

</div>

Thank you @BLI!

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 19, 2020, 9:13pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/4 "2020-04-19T21:13:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [April 20, 2020, 7:30am UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/5 "2020-04-20T07:30:04Z")

</div>

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](http://www.netlib.org/lapack/lug/node43.html).

---

<div class="post-metadata">

**Author:** ![Alex\_Ellison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alex_ellison/32/5749_2.png) [@Alex\_Ellison](https://discourse.julialang.org/u/Alex_Ellison)\
**Post date:** [April 20, 2020, 4:20pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/6 "2020-04-20T16:20:53Z")

</div>

Would you mind sharing how you’d do that?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [April 20, 2020, 5:12pm UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/7 "2020-04-20T17:12:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [April 21, 2020, 12:24am UTC](https://discourse.julialang.org/t/underdetermined-linear-system-general-solution/37865/8 "2020-04-21T00:24:54Z")

</div>

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

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