Numerical methods. Note 5

[~Home~]

Linear Least-Squares Problem

Numerical Methods. Note~ « 5 »

Least-Squares Solution. Let A be an n×m matrix, c 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|²

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 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 matrix A (with the help of QR decomposition).
Make a linear fit to the following data
x_i = i + ^sin(i)/₂, y_i = i + cos(i²), σ_i=0.1, i=0..10
Calculate the covariance matrix.