# LA - Bonus 2
[toc]

## Problem
Take *m = 50*, *n = 12*. Using MATLAB's `linspace`, define *t* to be the m-vector corresponding to linearly spaced grid points from 0 to 1.

Using MATLAB's `vander` and `fliplr`, define *A* to be the m x n matrix associated with least squares fitting on this grid by a polynomial of degree n - 1. Take b to be the function cos(4t) evaluated on the grid.

Now, calculate and print (to sixteen-digit precision) the least squares coefficient vector x by six methods:
1. Formation and solution of the normal equations, using MATLAB's `\`
2. QR factorization computed by mgs (modified Gram-Schmidt)
3. QR factorization computed by house (Householder triangularization)
4. QR factorization computed by MATLAB's qr (also Householder triangularization)
5. x = A \b in MATLAB (also based on QR factorization)
6. SVD, using MATLAB's svd.

The calculations above will produce six lists of twelve coefficients. In each list, shade with red pen the digits that appear to be wrong (affected by rounding error). Comment on what differences you observe. Do the normal equations exhibit instability?

## Code
### bonus2.m

% Define the Variables
m = 50;
n = 12;
t = linspace(0, 1, m);

% Define the Matrix A
A = fliplr(vander(t));
A = A(:, 1:n); % Take only the first n columns

% Define the Vector b
b = cos(4*t)';


% (a) Normal Equations
x_a = (A'*A) \ (A'*b);

% (b) QR Factorization using Modified Gram-Schmidt
[Q, R] = mgs(A);
x_b = R \ (Q' * b);

% (c) QR Factorization using Householder Triangularization
[Q, R] = house(A);
x_c = R \ (Q' * b);

% (d) QR Factorization using MATLAB's 'qr'
[Q, R] = qr(A, 0); % Economy size
x_d = R \ (Q' * b);

% (e) Using MATLAB’s Backslash Operator A\b
x_e = A \ b;

% (f) Singular Value Decomposition (SVD)
[U, S, V] = svd(A, 'econ');
x_f = V * (S \ (U' * b));

% Print the Results
format long g
disp('Normal Equations');
disp(x_a);
disp('QR (MGS)');
disp(x_b);
disp('QR (Householder)');
disp(x_c);
disp('QR (MATLAB)');
disp(x_d);
disp('MATLAB \\ operator');
disp(x_e);
disp('SVD');
disp(x_f);

### mgs.m

function [Q, R] = mgs(A)
    [m, n] = size(A);
    Q = zeros(m, n);
    R = zeros(n, n);
    for k = 1:n
        R(k, k) = norm(A(:, k));
        Q(:, k) = A(:, k) / R(k, k);
        for j = k+1:n
            R(k, j) = Q(:, k)' * A(:, j);
            A(:, j) = A(:, j) - Q(:, k) * R(k, j);
        end
    end
end

### house.m

function [Q, R] = house(A)
    [m, n] = size(A);
    Q = eye(m);
    R = A;
    for k = 1:n
        x = R(k:m, k);
        e = zeros(length(x), 1); e(1) = norm(x);
        u = x - e;
        u = u / norm(u);
        R(k:m, k:n) = R(k:m, k:n) - 2 * u * (u' * R(k:m, k:n));
        Q(k:m, :) = Q(k:m, :) - 2 * u * (u' * Q(k:m, :));
    end
    Q = Q';
end

## Result

### Normal Equations
1. 1.00000000==446473==
1. ==-1.73588168950921e-06==
1. -7.9999==2733282061==
1. -0.00==116408663073312==
1. 10.6==762208723219==
1. -0.0==456691644612114==
1. -5.==55385364161096==
1. -0.==250670185502446==
1. ==1.90558352040681==
1. -0.==153114645765417==
1. -0.3==06716135166133==
1. 0.0==756689052844102==

### QR (MGS)
1. 1.00000000==103769==
1. -4.2==6525959264673e-07==
1. -7.9999812==2547994==
1. -0.00031==7425649716468==
1. 10.6694==111405976==
1. -0.013==6954193525359==
1. -5.647==51845509882==
1. -0.07==43625403686041==
1. 1.69==233141281816==
1. 0.00==706880889388241==
1. -0.374==710650596597==
1. 0.088==1311602469381==

### QR (Householder)
1. 1.000000000==99661==
1. -4.2274==3233624982e-07==
1. -7.9999812356==7958==
1. -0.000318763==312221441==
1. 10.66943079==64629==
1. -0.0138202==903666128==
1. -5.6470756==2098654==
1. -0.0753160==354435472==
1. 1.693606976==14016==
1. 0.006032==09971610138==
1. -0.374241==699764291==
1. 0.08804057==54185232==

### QR (MATLAB)
1. 1.00000000099661
1. -4.2274287==7404515e-07==
1. -7.9999812356==8991==
1. -0.000318763==180767522==
1. 10.66943079==55492==
1. -0.0138202==865207034==
1. -5.6470756==3133516==
1. -0.0753160==172382638==
1. 1.6936069==5529251==
1. 0.006032==11469021733==
1. -0.374241==705889927==
1. 0.08804057==65072875==

### MATLAB `\` operator
1. 1.00000000099661
1. -4.2274==3364218306e-07==
1. -7.9999812356==7615==
1. -0.000318763==346323292==
1. 10.66943079==66411==
1. -0.0138202==909146187==
1. -5.6470756==1995939==
1. -0.07531603==65894195==
1. 1.693606976==80362==
1. 0.006032==09964510418==
1. -0.374241==699881279==
1. 0.08804057==54623562==

### SVD
1. 1.00000000099661
1. -4.2274287==4708903e-07==
1. -7.9999812356==9013==
1. -0.000318763==177247201==
1. 10.66943079==55211==
1. -0.0138202==863885616==
1. -5.6470756==3172688==
1. -0.0753160==164848571==
1. 1.6936069==5435502==
1. 0.006032==11541764992==
1. -0.374241==706209554==
1. 0.08804057==65679626==

## Observations
- **Normal Equations**\
==有很明顯的rounding errors和instability==，且位數越小越明顯。這種不穩定性是因為法方程涉及矩陣條件數的平方，放大了有限數值精度引起的誤差。
* **QR (MGS)**\
有些微的rounding errors但整體而言是比normal equations更穩定。
* **QR (Householder) and QR (MATLAB)**\
結果非常一致，幾乎沒有rounding errors。這裡產生的係數與 MATLAB 的內建 QR 和 SVD 幾乎相同，表示這兩個方法具有很強的數值穩定性。
- **MATLAB `\` operator**\
結果幾乎與 QR 和 SVD 一樣，穩定性算是很高。
- **SVD**\
奇異值分解法得到的結果與 QR 的非常接近，==提供了最穩定且精確的係數==，幾乎沒有rounding errors，是解決最小二乘問題的最佳選擇。