Least Squares and Pseudoinverse

Given an overdetermined system $\mathbf{A}\mathbf{x} \approx \mathbf{b}$ with $\mathbf{A} \in \mathbb{C}^{m \times n}$, $m > n$, the least-squares solution minimizes the residual:

$$\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x} - \mathbf{b}\|_2^2$$

The solution satisfies the normal equations: $$\mathbf{A}^H\mathbf{A}\hat{\mathbf{x}} = \mathbf{A}^H\mathbf{b}$$

x_hat, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)

lstsq returns:

• x_hat: the least-squares solution
• residuals: $\|\mathbf{A}\hat{\mathbf{x}} - \mathbf{b}\|_2^2$ (if $m > n$)
• rank: effective rank of $\mathbf{A}$
• sv: singular values of $\mathbf{A}$

The Moore-Penrose pseudoinverse $\mathbf{A}^\dagger$ generalizes the inverse to rectangular and rank-deficient matrices. If $\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^H$, then:

$$\mathbf{A}^\dagger = \mathbf{V}\boldsymbol{\Sigma}^\dagger\mathbf{U}^H$$

where $\Sigma^\dagger_{ii} = 1/\sigma_i$ for $\sigma_i > 0$ and $0$ otherwise.

A_pinv = np.linalg.pinv(A)
x_hat = A_pinv @ b   # equivalent to lstsq

For full-rank $\mathbf{A}$:

• $m > n$: $\mathbf{A}^\dagger = (\mathbf{A}^H\mathbf{A})^{-1}\mathbf{A}^H$ (left inverse)
• $m < n$: $\mathbf{A}^\dagger = \mathbf{A}^H(\mathbf{A}\mathbf{A}^H)^{-1}$ (right inverse)

When $\mathbf{A}$ is ill-conditioned, OLS amplifies noise. Tikhonov regularization adds a penalty:

$$\hat{\mathbf{x}}_{\mathrm{Tik}} = \arg\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x} - \mathbf{b}\|_2^2 + \alpha\|\mathbf{x}\|_2^2$$

The closed-form solution is: $$\hat{\mathbf{x}}_{\mathrm{Tik}} = (\mathbf{A}^H\mathbf{A} + \alpha\mathbf{I})^{-1}\mathbf{A}^H\mathbf{b}$$

alpha = 0.1
x_tik = np.linalg.solve(A.conj().T @ A + alpha * np.eye(n),
                         A.conj().T @ b)

This is exactly the MMSE estimator in wireless: with $\alpha = \sigma_n^2$, Tikhonov regularization gives the minimum mean-squared error (MMSE) estimate.

The QR factorization $\mathbf{A} = \mathbf{Q}\mathbf{R}$ where $\mathbf{Q}$ is unitary and $\mathbf{R}$ is upper triangular provides a numerically stable way to solve least squares:

$$\hat{\mathbf{x}} = \mathbf{R}^{-1}\mathbf{Q}^H\mathbf{b}$$

from scipy.linalg import qr, solve_triangular

Q, R = qr(A, mode='economic')   # reduced QR
x_qr = solve_triangular(R, Q.conj().T @ b)

QR is more numerically stable than the normal equations $(\mathbf{A}^H\mathbf{A})\mathbf{x} = \mathbf{A}^H\mathbf{b}$ because it avoids forming $\mathbf{A}^H\mathbf{A}$ (which squares the condition number).

When both $\mathbf{A}$ and $\mathbf{b}$ are noisy, ordinary least squares is biased. Total least squares minimizes:

$$\min_{\Delta\mathbf{A}, \Delta\mathbf{b}} \|[\Delta\mathbf{A} \mid \Delta\mathbf{b}]\|_F \quad \text{s.t.} \quad (\mathbf{A} + \Delta\mathbf{A})\mathbf{x} = \mathbf{b} + \Delta\mathbf{b}$$

The solution uses the SVD of $[\mathbf{A} \mid \mathbf{b}]$:

C = np.column_stack([A, b])
U, s, Vh = np.linalg.svd(C, full_matrices=False)
V = Vh.conj().T
x_tls = -V[:n, -1] / V[-1, -1]