# Systems of linear equations
The solution of a system of linear equations $Ax=b$ is obtained as


```python
import numpy as np

A = np.array([[1, 2], [3, -1]])
b = np.array([-1, 4])
```


```python
x = np.linalg.solve(A,b)
print(x)
```

    [ 1. -1.]


To check that this really is a solution, you should always compute the residual, which should vanish (to machine precision):


```python
r = A.dot(x)-b
print(r)
```

    [4.4408921e-16 0.0000000e+00]


The same solution is obtained using the command


```python
x = np.linalg.inv(A).dot(b)
print(x)
```

    [ 1. -1.]


The difference is that this requires computing the inverse of $A$, which is time consuming, whereas `solve(A,b)` applies a fast algorithm for solving linear equations, which does not require the inverse explicitly. This becomes relevant for large systems:


```python
import timeit
print(timeit.timeit(stmt='np.linalg.solve(A,b)',setup='import numpy as np; A = np.random.rand(1000,1000); b = np.random.rand(1000,1)',number=100))
print(timeit.timeit(stmt='np.linalg.inv(A).dot(b)',setup='import numpy as np; A = np.random.rand(1000,1000); b = np.random.rand(1000,1)',number=100))
```

    0.834245866
    2.5875389249999996


You can either use `linalg.inv()` and `linalg.dot()` methods in chain to solve a system of linear equations, or you can simply use the `solve()` method. The `solve()` method is the preferred way.

## Overdetermined systems
A system of linear equations is considered overdetermined if there are more equations than unknowns. For example, we have the overdetermined system

$$x_1 + x_2 =2 \\
x_1 = 1\\
x_2 = 0$$

In practice, we have a system $Ax=b$ where $A$ is a $m$ by $n$ matrix and $b$ is a $m$ dimensional vector, but $m$ is greater than $n$. In this case, the vector $b$ cannot be expressed as a linear combination of the columns of $A$. Hence, we can't find $x$ so that satisfies the problem $Ax=b$ (except in specific cases) but it is possible to determine $x$ so that $Ax$ is as close to $b$ as possible. So we wish to find $x$ which minimizes  $\begin{Vmatrix}Ax-b\end{Vmatrix}$. Considering the [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of $A$, we have that $Ax=b$ becomes $QRx=b$. Multiplying by $Q^T$ we obtain $Q^TQRx=Q^Tb$, and since $Q^T$ is orthogonal (this means that Q^T*Q=I) we have $Rx=Q^Tb$.

Now, this is a well defined system, $R$ is an upper triangular matrix and $Q^T*b$ is a vector. More precisely $b$ is the orthogonal projection of $b$ onto the range of $A$ and $\begin{Vmatrix}Ax-b\end{Vmatrix}=\begin{Vmatrix}Rx-Q^Tb\end{Vmatrix}$.

The function `linalg.lstsq()` provided by numpy returns the least-squares solution to a linear system equation and is able to solve overdetermined systems. Let's compare the solutions of `linalg.lstsq()` with the ones computed using the QR decomposition:


```python
A = np.array([[1, 1], [1, 0],[0, 1]])
b = np.array([2, 1, 0])
x = np.linalg.lstsq(A,b,rcond=None)[0]
print(x)
```

    [1.33333333 0.33333333]



```python
Q,R = np.linalg.qr(A) # qr decomposition of A
Qb = np.dot(Q.T,b) # computing Q^T*b (project b onto the range of A)
x = np.linalg.solve(R,Qb) # solving R*x = Q^T*b
print(x)
```

    [1.33333333 0.33333333]


As we can see, the solutions are the same.

This is the vector for which the norm of the residual $\begin{Vmatrix}r\end{Vmatrix}$ becomes minimal:


```python
r =  A.dot(x)-b
print(r)
np.linalg.norm(r)
```

    [-0.33333333  0.33333333  0.33333333]





    0.5773502691896257



An even simpler example $x=0$, $x=1$. Here the norm of the residual $\begin{Vmatrix}Ax-b\end{Vmatrix}$ is

$$\begin{Vmatrix}r\end{Vmatrix}=\sqrt{x^2+(x-1)^2}$$

and minimising this function (by finding the zero of the derivative) yields $x=1/2$. This is indeed what numpy returns:


```python
A = np.array([[1], [1]]);
b = np.array([0, 1]);
x = np.linalg.lstsq(A,b,rcond=None)[0]
print(x)
```

    [0.5]


## Underdetermined systems
As an example consider
$$x_1 + 2 x_2 + 3 x_3 + 4 x_4  = 1\\
5 x_1 + 6 x_2 + 7 x_3 + 8 x_4 = 2$$
This has an infinite number of solutions. `lstsq` returns one solution:


```python
A = np.array([[1,2,3,4],[5,6,7,8]])
b = np.array([1,2])
x = np.linalg.lstsq(A,b,rcond=None)[0]
print(x)
```

    [-0.05   0.025  0.1    0.175]


Using `scipy.optimize.nnls`, a non-negative least squares solver. Solve $argmin_x \begin{Vmatrix}Ax - b\end{Vmatrix}_2$ for $x\ge0$, returns a solution with as many components as possible equal to zero. 


```python
from scipy.optimize import nnls 
x, rnorm = nnls(A,b)
print(x)
print(rnorm)
```

    [1.65502277e-16 0.00000000e+00 0.00000000e+00 2.50000000e-01]
    0.0


To find all solutions, we need to determine the kernel of $A$. The nullspace or kernel of a matrix $A$ (denoted $\ker A$) is the set of all vectors $x$, for which $Ax=0$. If $x$ and $y$ are in the nullspace, then $c_1x+c_2y$ is also in the nullspace as 

$$A(c_1x+c_2y)=c_1(Ax)+c_2(Ay)=0+0=0$$

The nullspace is a vector space. When $A$ is viewed as a linear transformation, the nullspace is the subspace of $\mathcal{R}^n$ that is sent to 0 under the map $A$, hence the term "fundamental subspace."

An orthonormal basis $N=(n_1,\dots,n_k)$ of the kernel is returned by


```python
from scipy.linalg import null_space, orth
N = null_space(A)
print(N)
```

    [[-0.40008743 -0.37407225]
     [ 0.25463292  0.79697056]
     [ 0.69099646 -0.47172438]
     [-0.54554195  0.04882607]]


The nullspace consists of all vectors $x$ such that $Ax=0$. This defines a system of linear equations that can be solved to give the family of solutions

$$$$


```python
n0=N[:,0]
n1=N[:,1]
assert (abs(A.dot(n0))<=1e-14).all(), "Ax=0 for all x in nullspace {}".format(A.dot(n1))
assert (abs(A.dot(n1))<=1e-14).all(), "Ax=0 for all x in nullspace {}".format(A.dot(n1))
```

which defines a vector space with basis $\{n_0,n_1\}$. As there are two vectors in this basis, the dimension of the nullspace is 2.


```python
# example fundamental subspaces
A = np.array([[1,2,3,3],[2,0,6,2],[3,4,9,7]])
# column/row space (image)
print('column space: {}'.format(orth(A)))
print('row space: {}'.format(orth(A.T)))
# nullspace/ left nullspace (kernel)
print('nullspace: {}'.format(null_space(A)))
print('left nullspace: {}'.format(null_space(A.T)))
assert (abs(A.dot(null_space(A)))<=1e-14).all(), "Ax=0 for all x in nullspace {}".format(A.dot(null_space(A)))
# Fundamental Theorem of Linear Algebra: rank-nullity theorem (relates the dimensions of the four fundamental subspaces)
assert np.linalg.matrix_rank(A)==len(orth(A)[0])==len(orth(A.T)[0]), "The column and row spaces of an m×n matrix A both have dimension r, the rank of the matrix."
assert A.shape[1]-np.linalg.matrix_rank(A)==len(null_space(A)[1]), "The nullspace has dimension n−r."
assert A.shape[0]-np.linalg.matrix_rank(A.T)==len(null_space(A.T)[1]), "The left nullspace has dimension m−r."
# Fundamental Theorem of Linear Algebra: orthogonal spaces (the dot product v⋅w is 0)
assert (abs(null_space(A).T.dot(orth(A.T)))<=1e-14).all(), "The nullspace and row space are orthogonal."
assert (abs(null_space(A.T).T.dot(orth(A)))<=1e-14).all(), "The left nullspace and the column space are also orthogonal."
# Fundamental Theorem of Linear Algebra:  orthonormal basis (singular value decomposition)
[U,s,V]=np.linalg.svd(A)
S = np.zeros(A.shape, dtype=complex)
S[:(A.shape[0]), :(A.shape[0])] = np.diag(s)
assert np.allclose(A, np.dot(U, np.dot(S, V))), "any matrix M can be written as dot product of an m x m unitary martix, an m x n matrix with nonnegative values in diagonal, and an b x n unitary matrix."
```

    column space: [[-0.31994238 -0.36841839]
     [-0.41936816  0.88119879]
     [-0.84956884 -0.29623738]]
    row space: [[-0.25358142  0.17588223]
     [-0.27620609 -0.66896915]
     [-0.76074427  0.52764668]
     [-0.52978752 -0.49308692]]
    nullspace: [[-0.94415867 -0.11543961]
     [ 0.03383545 -0.68923555]
     [ 0.32599804 -0.19126531]
     [-0.03383545  0.68923555]]
    left nullspace: [[ 0.87287156]
     [ 0.21821789]
     [-0.43643578]]


Using sympy to solve the equation set symbolically


```python
from sympy import * 

x_1, x_2, x_3, x_4 = symbols('x_1 x_2 x_3 x_4')

res = solve([Eq(1*x_1+2*x_2+3*x_3+4*x_4, 1),
             Eq(5*x_1+6*x_2+7*x_3+8*x_4, 2)],
             [x_1, x_2, x_3, x_4])
print(res)
```

    {x_1: x_3 + 2*x_4 - 1/2, x_2: -2*x_3 - 3*x_4 + 3/4}



```python

```

02613a47861c1d9d5f54ddb5f5fdb5e28588f925
ipynb
Jupyter Notebook
notebooks/Systems of linear equations.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
2022-01-08T23:52:57.000Z
2022-01-08T23:52:57.000Z
notebooks/Systems of linear equations.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
notebooks/Systems of linear equations.ipynb
OleBo/MathSo
1f9fa0492d467c0bb0479768c503eee7723ae777
MIT
2022-01-02T21:14:30.000Z
2022-01-02T21:14:30.000Z
Qwen/Qwen-72B
1. YES
2. YES