Octave for Chapter 9
Examples and Templates in this section provide sample Octave code for least squares, and -factorization. You can access our code through the link at the bottom of each template. Feel free to modify the code and experiment to learn more!
You can write your own code using Octave software or online Octave cells. To access Octave cells online, go to the Sage Math Cell Webpage, select OCTAVE as the language, enter your code, and press EVALUATE.
To ”save” or share your online code, click on the Share button, select Permalink, then copy the address directly from the browser window. You can store this link to access your work later or share this link with others. You will need to get a new Permalink every time you modify the code.
Octave Tutorial
Least Squares
The following theorem establishes one way to implement least squares.
th:bestApprox Let be an matrix, let be a column vector in . Consider the matrix equation
- (a)
- Any solution to the normal equations is a best approximation to a solution to in the sense that is minimized.
- (b)
- If the columns of are linearly independent, then is invertible and is given uniquely by
We can implement this process as follows.
We show how to find a least squares solution to .
% Define matrix A
A=[9 -3 3;
1 -1 1;
4 0 1;
1 1 1;
9 -1 2];
% Define vector b
b=[3;1;1;2;4];
% Now we solve for z
z=inv(transpose(A)*A)*transpose(A)*b
To do least squares in Octave, we use the backslash operator (). Depending on the type of matrix, backslash uses different techniques to compute the least squares solution more efficiently. For more information about inv and the backslash operator () see Reference.
The following two examples illustrate the two implementations of least squares side-by-side.
We will compare the results of the two implementations of least squares for a
coefficient matrix.
% Define matrix A
A=[4 -1 0;
2 -1 1;
6 3 1;
-1 1 2;
3 1 -2;
-2 8 7;
1 1 -5];
% Define vector b
b=[-2;4;1;2;-5;1;-3];
% Implement least squares using our original method
z1=inv(transpose(A)*A)*transpose(A)*b
% Implement least squares using the backslash operator
z2=A\b
We will compare the performance of two the implementations of least squares for a
coefficient matrix.
% Define matrix A
A=rand(2000,20); % rand(m,n) creates an m by n matrix whose entries are between 0 and 1
% Define vector b
b=rand(2000,1);
% Start the timer
tic
% Implement least squares using our original method
z1=inv(transpose(A)*A)*transpose(A)*b;
% Stop the timer
timeINV=toc
% Start the timer
tic
% Implement least squares using the backslash operator
z2=A\b;
% Stop the timer
timeBSlash=toc
Each time you run the algorithm, you will get slightly different run times for each implementation. Most of the time, you will see that the backslash operator outperforms our initial method.
The backslash operator () can be used to solve systems of equations when a solution exists. Here is an example.
Consider the system \begin{equation} \begin{array}{ccccccccc} x &- &y&&&&&= &0 \\ 2x& -&2y&+&z&+&2w&=&4\\ & &y&&&+&w&=&0\\ & &&&2z&+&w&=&5 \end{array} \end{equation}
(You may remember this system from Augmented Matrix Notation and Elementary Row Operations.) We will solve this system using the backslash operator ().
% Define the coefficient matrix A
A = [1 -1 0 0;
2 -2 1 2;
0 1 0 1;
0 0 2 1];
% Define vector b
b = [0;4;0;5];
% Solve the system
x = A \ b
% Compare the solution above to the solution obtained
% using rref
A_b=[A b];
rref(A_b)
While the backslash () operator can be used to solve feasible systems, we urge you to
be very cautious because
- If a system has infinitely many solutions, the backslash () operator will find one particular solution.
- If a system has no solutions, the backslash () operator will find an approximate solution using least-squares.
The way the answer is presented, you cannot tell the difference!
-factorization
Recall the definition of -factorization.
def:QR-factorization Let be an matrix with independent columns.
A QR-factorization of expresses it as where is with orthonormal columns
and is an invertible and upper triangular matrix with positive diagonal
entries.
Find the -factorization of
and demonstrate that matrix is orthogonal (columns are orthonormal).
% Define A
A=[2 -1 4 0;
3 1 1 -2;
1 5 -1 6;
-3 4 0 7];
% Find the QR-factorization of A
[Q,R]=qr(A)
% Verify that columns of Q are orthonormal
% You may get numbers that are close to zero but not zero
% due to round-off error
transpose(Q)*Q
% why does this verify that the columns are orthonormal?
In the code above, we used the product to verify that the columns of are
orthonormal. What do you expect this product to look like if the columns are
orthonormal? What will the product look like if they are not? You will be asked to
explain your reasoning in Problem prob_oct_orth.
Recall the following algorithm for using -factorization to approximate eigenvalues from -Factorization.
alg:qrEig Let be an invertible matrix.
Step 1: Define and factor it as .
Step 2: Define and factor it as .
Step 3: Define and factor it as .
Step : Define (where ).
Note that is similar to (in fact, ), and hence each has the same eigenvalues as . If the eigenvalues of are real and have distinct absolute values, the sequence of matrices converges to an upper triangular matrix with these eigenvalues on the main diagonal.
Compare the following example to Example QR-algortihm-2x2-025425.
Use the algorithm to approximate the eigenvalues of .
Here is a very simple implementation of the algorithm.
% Define A
A=[1 1;
2 0];
% Set the desired number of iterations
iter=20;
% Implement the algorithm. Note that we are over-writing the previous A with every iteration.
for i=1:iter
[Q,R]=qr(A);
A=R*Q
end
Observe that the main diagonal entries are approaching the eigenvalues of .
Octave Exercises
Attempt to solve each of the following systems of equations using the backslash ()
operator, and the rref function. Interpret your results.
- (a)
- \begin{equation} \begin{array}{ccccccccc} 3x &- &y&+&4z&= &2 \\ &&y&+&z&=&1\\ -2x&&&+&3z&=&1 \end{array} \end{equation}
What is true about this system?
The system is consistent. Using the backslash operator produces the same result as using the rref function. The system is inconsistent. The backslash operator gives the least squares approximation. The system is consistent but has infinitely many solutions. The backslash operator gives one particular solution. - (b)
- \begin{equation} \begin{array}{ccccccccc} 3x &- &y&+&4z&= &2 \\ x&-&y&+&2z&=&4\\ -2x&+&y&-&3z&=&-3 \end{array} \end{equation}
What is true about this system?
The system is consistent. Using the backslash operator produces the same result as using the rref function. The system is inconsistent. The backslash operator gives the least squares approximation. The system is consistent but has infinitely many solutions. The backslash operator gives one particular solution. - (c)
- \begin{equation} \begin{array}{ccccccccc} 3x &- &y&+&4z&= &2 \\ x&-&y&+&2z&=&1\\ -2x&+&y&-&3z&=&1 \end{array} \end{equation}
What is true about this system?
The system is consistent. Using the backslash operator produces the same result as using the rref function. The system is inconsistent. The backslash operator gives the least squares approximation. The system is consistent but has infinitely many solutions. The backslash operator gives one particular solution.
Use least squares to find a line of best fit for the points shown in the GeoGebra
interactive below. Enter the equation of the line into the GeoGebra interactive to
view your result.
- (a)
- Find the residual, , for your line of best fit.
- (b)
- The original line, initially shown in the interactive, visually appeared to be a pretty good fit. Find the residual for the line and compare it to the residual you got for your line. Is your line a better fit?
Explain how to interpret the product to determine whether matrix has
orthonormal columns. Prove your claims.
Modify the code in Example ex:qrEig to approximate the eigenvalues of . Compare your
answers to the answers you got for Problem prob:3x3fromKuttler1.
Modify the code in Example ex:qrEig for the matrix . What are you observing? Are we
getting close to finding the eigenvalues of ? Explain what is happening and
why.
2024-10-01 02:16:33