Showing that the pseudoinverse solution minimizes the sums of squares of the estimated parameters $B$ for $Y = X B + \epsilon$, in the case of a particular rank deficient matrix, from a tutorial by Cyril Pernet.


```
import numpy as np
import numpy.linalg as npl
```


```
# Setting up the design
a_block = np.array([-0.1, 0, 0.1])
baseline = 10
activation = 11
n_on_off = 3
# The data
on_off = np.hstack((a_block + baseline, a_block + activation))
Y = np.tile(on_off, (n_on_off,))
Y = Y[:, None] # Make Y into a column vector
# Regressors in the redundant design
x_on = np.hstack((np.zeros(len(a_block)), np.ones(len(a_block))))
x_off = 1 - x_on
X_over_part = np.column_stack((x_on, x_off, np.ones_like(x_on)))
X_over = np.tile(X_over_part, (n_on_off, 1))
Y, X_over
```




    (array([[  9.9],
            [ 10. ],
            [ 10.1],
            [ 10.9],
            [ 11. ],
            [ 11.1],
            [  9.9],
            [ 10. ],
            [ 10.1],
            [ 10.9],
            [ 11. ],
            [ 11.1],
            [  9.9],
            [ 10. ],
            [ 10.1],
            [ 10.9],
            [ 11. ],
            [ 11.1]]), array([[ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.],
            [ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.],
            [ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 0.,  1.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.],
            [ 1.,  0.,  1.]]))



The design matrix `X_over` is not full rank:


```
npl.matrix_rank(X_over)
```




    2



That means it has a *null space* - a set of vectors $k$ for which $X k = 0$


```
import sympy
sympy.init_printing(use_latex=True)
```


```
X = sympy.Matrix(X_over)
null_space = X.nullspace()
null_space
```

The null space is a list of all the basis vectors spanning the null space:


```
type(null_space)
```




    list



In our case, the matrix is only rank deficient by one column, and there is only one basis vector for the null space:


```
len(null_space)
```


```
null_vector = null_space[0]
null_vector
```

The null space is all vectors $k$ formed by multiplying a scalar $a$ times the null space basis vector `null_vector`.

This in turn means that if we find some parameter vector $B$ as the solution for $Y = X B$, then $Y = X B + X k = X (B + k)$ also.

For example, the parameters that went into forming $Y$ in our case are $1, 0, 10$ (activation, baseline, mean signal).  Not surprisingly, this vector minimizes the sum of squares for $\epsilon = Y - X B$:


```
beta_in = np.array([[1, 0, 10]]).T
Y - X_over.dot(beta_in)
```




    array([[-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1]])



But we get exactly the same fit by adding any multiple of `null_vector` to $B$:


```
null_numeric = np.array(null_vector).astype(float)
Y - X_over.dot(beta_in + null_numeric * 42)
```




    array([[-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1],
           [-0.1],
           [ 0. ],
           [ 0.1]])



This means that the residuals $\epsilon$ cannot distinguish between any of the possible parameter vectors $B + k$, where $k$ is a vector in the null space of $X$.

The pseudoinverse operator chooses one of the infinite possible solutions from $B + k$:


```
beta_hat = npl.pinv(X_over).dot(Y)
beta_hat
```




    array([[ 4.],
           [ 3.],
           [ 7.]])



In fact, this is the solution that minimizes the sum of squares of the chosen parameters in `beta_hat`:


```
# The fitted parameters can be beta_in (the input paramters) + a * null_vector
a = sympy.symbols('a')
fitted_params = beta_in + a * null_vector
fitted_params
```


```
# pinv minimized the sums of squares of the fitted parameters
ss = sympy.simplify(sum(e ** 2 for e in fitted_params))
ss
```


```
# What multiple `a` of the null vector minimizes the sums of squares?
a_min_ss = sympy.solve(sympy.diff(ss), a)[0]
a_min_ss
```


```
# That's what we get from the pseudoinverse model fit
sympy.Matrix(beta_in + a_min_ss * null_vector)
```

579a997d325bb6698d2f2069b179f3e0e70af69f
ipynb
Jupyter Notebook
pinv_minimizing.ipynb
matthew-brett/sketch-books
ed5afe08f2787880e70e30b023bc9fdf49338e3f
CC-BY-3.0
pinv_minimizing.ipynb
matthew-brett/sketch-books
ed5afe08f2787880e70e30b023bc9fdf49338e3f
CC-BY-3.0
pinv_minimizing.ipynb
matthew-brett/sketch-books
ed5afe08f2787880e70e30b023bc9fdf49338e3f
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES