Open in the MATLAB editor: f08bh_example
function f08bh_example fprintf('f08bh example results\n\n'); % Going to solve Least squares problem Ax = b, m>n. m = 6; n = 5; a = [-0.09, 0.14, -0.46, 0.68, 1.29; -1.56, 0.2, 0.29, 1.09, 0.51; -1.48, -0.43, 0.89, -0.71, -0.96; -1.09, 0.84, 0.77, 2.11, -1.27; 0.08, 0.55, -1.13, 0.14, 1.74; -1.59, -0.72, 1.06, 1.24, 0.34]; b = [ 7.4, 2.7; 4.2, -3.0; -8.3, -9.6; 1.8, 1.1; 8.6, 4.0; 2.1, -5.7]; % QR factorization of A with column pivoting = Q*(R1 R2 )*(P^T) % (0 R22) [qr, jpvt, tau, info] = f08bf( ... a, zeros(n,1,'int64')); % QRP'X = B, => RP'X = Q'B = C % Compute C = Q'B [c, info] = f08ag( ... 'Left', 'Transpose', qr, tau, b); % Determine the rank, k, of R relative to tol; % Choose tol to reflect the relative accuracy of the input data tol = 0.01; k = find(abs(diag(qr)) <= tol*abs(qr(1,1))); if numel(k) == 0 k = numel(diag(qr)); else k = k(1)-1; end fprintf('Tolerance used to estimate the rank of A\n %11.2e\n', tol); fprintf('Estimated rank of A\n %d\n', k); % Compute the RZ (TZ) factorization of the first k rows of (R1 R2) [rz, taurz, info] = f08bh( ... qr(1:k,:)); % Now, (TZ)P'X = C on first k rows of C % Let ZP'X = T^{-1}C = Y (on first k rows) y = zeros(n, 2); y(1:k, :) = inv(triu(rz(1:k,1:k)))*c(1:k,:); % ZP'X = Y => P'X = Z^T Y = W; Form W = Z^TY. [w, info] = f08bk( ... 'Left', 'Transpose', int64(n-k), rz, taurz, y); % P'X = W => X = PW, ' x = zeros(n, 2); for i=1:n x(jpvt(i), :) = w(i, :); end fprintf('\nLeast-squares solution(s)\n'); disp(x); % Compute estimates of the square roots of the residual sums of squares rnorm = [norm(c(k+1:6,1)), norm(c(k+1:6,2))]; fprintf('Square root(s) of the residual sum(s) of squares\n'); disp(rnorm);
f08bh example results Tolerance used to estimate the rank of A 1.00e-02 Estimated rank of A 4 Least-squares solution(s) 0.6344 3.6258 0.9699 1.8284 -1.4402 -1.6416 3.3678 2.4307 3.3992 0.2818 Square root(s) of the residual sum(s) of squares 0.0254 0.0365