Section B.7 Best approximation
Warning B.7.1.
This page contains several redefinitions of a Python function inprod
. You will get incorrect results or errors if you Evaluate
a Sage cell defining that function in one subsection below, and then Evaluate
Sage cells that use a function by that same name in a different subsection below without evaluating the appropriate Sage cell near the beginning of that different subsection that redefines that inprod
function to the appropriate formula used in the example in that particular subsection.
The best way to avoid this problem is to Evaluate
every Sage cell in a particular subsection, from the beginning, in order.
Subsection B.7.1 Approximating a matrix
Here we will use Sage to carry out the calculations in Example 38.3.1. We carry out the whole process, including Gram-Schmidt to obtain the orthogonal basis. The subspace U is described as consisting of [those] upper-triangular matrices with (1,2) entry equal to the trace of the matrix. Parametrically, we can describe U as
\begin{equation*}
U = \left\{ \begin{bmatrix} a \amp a + b \\ 0 \amp b \end{bmatrix} \right\} \text{,}
\end{equation*}
where a and b are free parameters. This leads to (non-orthogonal) basis
\begin{equation*}
\left\{
\begin{bmatrix} 1 \amp 1 \\ 0 \amp 0 \end{bmatrix},
\begin{bmatrix} 0 \amp 1 \\ 0 \amp 1 \end{bmatrix}
\right\} \text{,}
\end{equation*}
where these matrices correspond to parameter choices \{a = 1, b = 0\} and \{a = 0, b = 1\}\text{,} respectively.
Set up.
First let's load our initial basis vectors for U into Sage.xxxxxxxxxx
A1 = matrix( [ [1, 1], [0, 0] ] )
A2 = matrix( [ [0, 1], [0, 1] ] )
print(A1)
print()
print(A2)
xxxxxxxxxx
def inprod(A,B):
return (B.T * A).trace()
print("inprod function loaded")
B.T
is a Sage βshortcutβ for B.transpose()
.Gram-Schmidt That Thing.
It's the dead duo's time to shine.xxxxxxxxxx
E1 = A1
E2 = A2 - (inprod(A2,E1) / inprod(E1,E1)) * E1
print(E1)
print()
print(E2)
xxxxxxxxxx
E2 = 2 * E2
print(E2)
Projection time.
The point of this example is to compute the matrix in U that is closest to being the identity matrix. For that, we want to compute \proj_U I\text{.}xxxxxxxxxx
ID = identity_matrix(2)
PROJ = (inprod(ID,E1) / inprod(E1,E1)) * E1 + (inprod(ID,E2) / inprod(E2,E2)) * E2
print(PROJ)
xxxxxxxxxx
O_PROJ = ID - PROJ
inprod(PROJ, O_PROJ)
xxxxxxxxxx
sqrt(inprod(O_PROJ, O_PROJ))
Subsection B.7.2 Approximating a function
Here we will use Sage to carry out the calculations in Example 38.3.2. As we have already demonstrated using Sage to apply Gram-Schmidt on our initial basis for this problem in Section B.6, we'll skip that part this time and proceed immediately to entering in our function, our orthogonal basis, and our inner product function.xxxxxxxxxx
f = sin(pi * x)
print(f)
e1 = 1
e2 = x - 1/2
e3 = x^2 - x + 1/6
print(e1,',',e2,',',e3)
def inprod(f,g):
return integrate(f*g,x,0,1)
print("inprod function loaded")
xxxxxxxxxx
proj = \
(inprod(f,e1) / inprod(e1,e1)) * e1 \
+ (inprod(f,e2) / inprod(e2,e2)) * e2 \
+ (inprod(f,e3) / inprod(e3,e3)) * e3
print(proj)
xxxxxxxxxx
o_proj = f - proj
err = sqrt(inprod(o_proj,o_proj))
print(err.n())
xxxxxxxxxx
naive_err_vec = f - (4*x - 4*x^2)
naive_err = sqrt(inprod(naive_err_vec,naive_err_vec))
print(naive_err.n())