Numerical methods. Note 3

[ Home ]

Linear Least-Squares Problem

Numerical Methods. Note « 3 »

Least-Squares Solution. Let A be an n×m matrix, let c be an m-component vector and b an n-component vector. For n > m the equation Ac = b generally has no solution. However one can find the "best possible" (least squares) solution which minimizes the Eucledian norm of the difference between Ac and b. The problem can be solved by the QR decomposition. The matrix A factorizes as A=QR, where Q is n×m matrix with orthogonal columns and R is an m×m upper triangular matrix. We can then minimize

|Ac-b|² = |QRc-b|² = |Rc - Q^Tb|² + |(1 - QQ^T)b|² ≥ |(1 - QQ^T)b|²

by solving a m×m set of linear equations Rc = Q^Tb by simple back-substitution.

Linear Least-Squares Fit is a problem of fitting n data points y_i (with errors σ_i) by a linear combination of m functions F(x) = ∑_k c_k f_k(x). The best fit is then when the square deviation χ² = ∑_i (^{[y_i - F(x_i)]}/_{σ_i})² is minimized. One can recognise the above problem of minimizing |Ac - b|² where A_ik = ^f_k(x_i)/_{σ_i} , b_i = ^y_i/_{σ_i} with the formal solution c = R^-1Q^Tb (in practice you should backsubstitute Rc = Q^Tb).

The error of the least-squares fit Δc_i can be approximated as

Δc_i = ∑_k ^∂c_i/_{∂y_k}×σ_k = ∑_k ^∂c_i/_{∂b_k} ≈ √( ∑_k [^∂c_i/_{∂b_k}]² ).

The error matrix (also called the covariance matrix) E_ij = ⟨Δc_iΔc_j⟩ is then simply

E = (^∂c/_∂b) (^∂c/_∂b)^T = R^-1(R^-1)^T = A^-1(A^-1)^T = (R^TR)^-1 = (A^TA)^-1

The diagonal elements of the covarinace matrix are the squares of the errors in the corresponding coefficients Δc_i = √(E_ii). The off-diagonal elements characterise correlations in the data -- if the (normalised) off-diagonal element is close to one, ^E_ij/_{√(E_iiE_jj)} ≈ 1, then the coefficients c_i and c_j can not be reliably estimated from the data.

The inverse A^-1 of a matrix A can be calculated by solving m linear systems A_ikz_k^( j ) = δ_ij, j=1:m. The inverse matrix then will be equal (A^-1)_ij = z_i^( j ).

Problems

Make a subroutine that calculates an inverse of a given matrix A (with the help of QR decomposition).
Make a subroutine that fits a given data-set [x_i, y_i, σ_i, i=1..n] by a linear combination of functions [f_k(x), k=1..m]. The subroutine must return the vector of coefficients c and the covariance matrix.
Make a linear fit to the following data
x_i = i + ^sin(i)/₂, y_i = i + cos(i²), σ_i = sin(i+1)², i=0..10
(husk, at i = 0..10, ikke 1..10)