I don't know what I'm wrong..

years = [1000, 1650, 1800, 1900, 1950, 1960, 1970, 1980, 1990];
population = [0.34, 0.545, 0.907, 1.61, 2.51, 3.15, 3.65, 4.2, 5.3];

x = years;
y = population;
A = [ones(length(x), 1), x'];
b = y' ./ (1 + x');

[Q, R] = householderQR(A);

c1 = Q' * b;
c1 = c1(1:2);
x_hh = R \ c1;

tt = linspace(min(x), max(x), 100);
yy = (x_hh(1) + x_hh(2)*tt) ./ (1 + tt);

clf;
hold on;
plot(tt, yy, 'LineWidth', 2);
plot(x, y, 'k*', 'MarkerSize', 10, 'LineWidth', 2);
hold off;
grid on;
legend('Fitting Curve', 'Data Points');
title('World Population Fitting Using Householder QR Factorization');
xlabel('Year');
ylabel('Population (billions)');

function [Q, R] = householderQR(A)
[m, n] = size(A);
Q = eye(m);
R = A;
for k = 1:n
    % Step 1: Define alpha
    alpha = -sign(R(k, k)) * norm(R(k:m, k));
    
    % Step 2: Define v
    v = R(k:m, k);
    v(1) = v(1) - alpha;
    v = v / norm(v);
    
    % Step 3: Define H
    H = eye(m-k+1) - 2 * (v * v');
    H_full = eye(m);
    H_full(k:m, k:m) = H;
    
    % Step 4: Update R and Q
    R = H_full * R;
    Q = Q * H_full';
end
% Adjust R to be upper triangular of size n x n
R = R(1:n, :);
end

What I wrote is the code of the below question.

A quick way of finding a function of the form f(x) ≈ (a + bx) / (1 + cx) is to apply the least squares method to the problem (1+cx)f(x) ≈ (a+bx). Use this technique to fit the world population (billions) data given using Householder QR method (without using built-in functions). Make a plot of the original data points along with resulting curve.

I don't know why my code of my plot is not working. Could someone modify my code to work well?

[Sathvik]

What parts of the code were you asked to modify as part of the question?

[Sam Chak]

Hi @matlabnm

Based on the information from Wikipedia, the algorithm for Householder transformation suggests that the computation of alpha (α) should start from the 2nd element. However, in your for-loop, it appears to compute alpha from the 1st element. This deviation from the recommended procedure might affect the results.

Additionally, the procedure you are following involves a redefined signum function where is equal to 1. However, in MATLAB, the sign(0) function returns zero. To address this discrepancy, you may need to artificially set the value of the signum function sign(0) at the origin to ensure consistency with the desired algorithm.

sign(0)

ans = 0

function [Q, R] = householderQR(A)
    [m, n]  = size(A);
    Q       = eye(m);
    R       = A;
    
    for k = 1:n
        % Step 1: Define alpha
        alpha   = -sign(R(k, k)) * norm(R(k:m, k));
    
        % Step 2: Define v
        v       = R(k:m, k);
        v(1)    = v(1) - alpha;
        v       = v / norm(v);
    
        % Step 3: Define H
        H       = eye(m-k+1) - 2 * (v * v');
        H_full  = eye(m);
        H_full(k:m, k:m) = H;
    
        % Step 4: Update R and Q
        R       = H_full * R;
        Q       = Q * H_full';
    end
    
    % Adjust R to be upper triangular of size n x n
    R       = R(1:n, :);
end