Je bent je ingevulde velden bij deze pagina aan het verwijderen. Ben je zeker dat je dit wilt doen?
You are erasing your filled-in fields on this page. Are you sure that is what you want?
Nieuwe Versie BeschikbaarNew Version Available
Er is een update van deze pagina. Als je update naar de meest recente versie, verlies je mogelijk je huidige antwoorden voor deze pagina. Hoe wil je verdergaan ?
There is an updated version of this page. If you update to the most recent version, then your current progress on this page will be erased. Regardless, your record of completion will remain. How would you like to proceed?
Examples and Templates in this section provide sample Octave code for least squares, and \(QR\)-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.
?? Let \(A\) be an \(m\times n\) matrix, let \(\vec {b}\) be a column vector in \(\RR ^m\). Consider the matrix equation
\[A\vec {x}=\vec {b}\]
1.
Any solution \(\vec {z}\) to the normal equations
\[\left (A^TA\right )\vec {z}=A^T\vec {b}\]
is a best approximation to a solution to \(A\vec {x}=\vec {b}\) in the sense that \(\norm {\vec {b}-A\vec {z}}\) is minimized.
2.
If the columns of \(A\) are linearly independent, then \(A^TA\) is invertible and \(\vec {z}\) is given uniquely by
\[\vec {z}=\left (A^TA\right )^{-1}A^T\vec {b}\]
We can implement this process as follows.
We show how to find a least squares solution to \(A\vec {x}=\vec {b}\).
% 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
Finding matrix inverses is both computationally expensive and unstable.
To do least squares in Octave, we use the backslash operator (\(\setminus \)). 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
(\(\setminus \)) 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 \(7\times 3\) 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 \(2000\times 20\) 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 (\(\setminus \)) can be used to solve systems of equations when a solution exists. Here is an example.
% 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 (\(\setminus \)) 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 (\(\setminus \)) operator will find one particular solution.
If a system has no solutions, the backslash (\(\setminus \)) operator will find an approximate solution using least-squares.
The way the answer is presented, you cannot tell the difference!
\(QR\)-factorization
Recall the definition of \(QR\)-factorization.
?? Let \(A\) be an \(m \times n\) matrix with independent columns. A QR-factorization of \(A\) expresses it as \(A = QR\)
where \(Q\) is \(m \times n\) with orthonormal columns and \(R\) is an invertible and upper triangular matrix with positive diagonal entries.
and demonstrate that matrix \(Q\) 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 \(Q^TQ\) to verify that the columns of \(Q\) 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 .
Recall the following algorithm for using \(QR\)-factorization to approximate eigenvalues from \(QR\)-Factorization.
?? Let \(A\) be an invertible matrix.
Step 1: Define \(A_{1} = A\) and factor it as \(A_{1} = Q_{1}R_{1}\).
Step 2: Define \(A_{2} = R_{1}Q_{1}\) and factor it as \(A_{2} = Q_{2}R_{2}\).
Step 3: Define \(A_{3} = R_{2}Q_{2}\) and factor it as \(A_{3} = Q_{3}R_{3}\).
Note that \(A_{k + 1}\) is similar to \(A_{k}\) (in fact, \(A_{k+1} = R_{k}Q_{k} = (Q_{k}^{-1}A_{k})Q_{k}\)), and hence each \(A_{k}\) has the same eigenvalues as \(A\). If the eigenvalues of \(A\) are real and have
distinct absolute values, the sequence of matrices \(A_{1}, A_{2}, A_{3}, \dots \) converges to an upper triangular matrix with these eigenvalues on the
main diagonal.
Compare the following example to Example ??.
Use the \(QR\) algorithm to approximate the eigenvalues of \(A=\begin{bmatrix}1 & 1\\2 & 0\end{bmatrix}\).
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 \(A\).
Octave Exercises
Attempt to solve each of the following systems of equations using the backslash (\(\setminus \)) operator, and the rref function. Interpret
your results.
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.
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.
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.
1.
Find the residual, \(\norm {\vec {b}-A\vec {z}}\), for your line of best fit.
2.
The original line, \(y=x+1\) initially shown in the interactive, visually appeared to be a pretty good fit. Find the residual
for the line \(y=x+1\) and compare it to the residual you got for your line. Is your line a better fit?
Explain how to interpret the product \(Q^TQ\) to determine whether matrix \(Q\) has orthonormal columns. Prove your claims.
Modify the code in Example 12 to approximate the eigenvalues of \(A=\begin{bmatrix} 9 & 2 & 8\\ 2 & -6 & -2\\ -8 & 2 & -5\end{bmatrix}\). Compare your answers to the answers you got for
Problem ??.
Modify the code in Example 12 for the matrix \(A=\begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix}\). What are you observing? Are we getting close to finding the eigenvalues of \(A\)?
Explain what is happening and why.
Linear Algebra: An Interactive Introduction