# Chapter 6 - Orthogonality

## Section 6.1

There are many ways for things to be orthogonal.

vector([1,2,3]).dot_product(vector([-3,0,1]))
 0 0
vector([1,2,3])*vector([-3,0,1])
 0 0

For instance, we have some polynomials that are, if you consider the inner product to be $\int_{-1}^1 f(x)g(x)dx$.

for i in range(4): show(expand(legendre_P(i,x)))
for i in range(4): for j in range(i+1,4): f(x) = legendre_P(i,x) g(x) = legendre_P(j,x) html("$$\int_{-1}^1 \left(%s\\right) \left(%s\\right)\; dx = %s$$"%(latex(f(x)),latex(g(x)),integrate(legendre_P(i,x)*legendre_P(j,x),(x,-1,1))))
 \int_{-1}^1 \left(1\right) \left(x\right)\; dx = 0 \int_{-1}^1 \left(1\right) \left(\frac{3}{2} \, {\left(x - 1\right)}^{2} + 3 \, x - 2\right)\; dx = 0 \int_{-1}^1 \left(1\right) \left(\frac{5}{2} \, {\left(x - 1\right)}^{3} + \frac{15}{2} \, {\left(x - 1\right)}^{2} + 6 \, x - 5\right)\; dx = 0 \int_{-1}^1 \left(x\right) \left(\frac{3}{2} \, {\left(x - 1\right)}^{2} + 3 \, x - 2\right)\; dx = 0 \int_{-1}^1 \left(x\right) \left(\frac{5}{2} \, {\left(x - 1\right)}^{3} + \frac{15}{2} \, {\left(x - 1\right)}^{2} + 6 \, x - 5\right)\; dx = 0 \int_{-1}^1 \left(\frac{3}{2} \, {\left(x - 1\right)}^{2} + 3 \, x - 2\right) \left(\frac{5}{2} \, {\left(x - 1\right)}^{3} + \frac{15}{2} \, {\left(x - 1\right)}^{2} + 6 \, x - 5\right)\; dx = 0  

Or you can try some complex vectors.

vector([1+i,2])*vector([1-i,4])
 0 0

The functions $\sin(kx)$ and $\cos(\ell x)$ are also all orthogonal on $[0,2\pi]$.

@interact def _(i=([1..5]),j=([1..5])): html("The integral $$\int_0^{2\pi} \sin(%s x)\sin(%s x)\; dx$$"%(i,j)) html("is $%s$"%latex(integrate(sin(i*x)*sin(j*x),(x,0,2*pi)))) html("The integral $$\int_0^{2\pi} \cos(%s x)\sin(%s x)\; dx$$"%(i,j)) html("is $%s$"%latex(integrate(cos(i*x)*sin(j*x),(x,0,2*pi))))

## Click to the left again to hide and once more to show the dynamic interactive window

Anyway, now we have a way to talk about what we claimed before - that the row space and null space are "perpendicular"!

M = matrix([[1,2,3],[1,2,3],[1,2,3]]); show(M)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 2 & 3 \\ 1 & 2 & 3 \\ 1 & 2 & 3 \end{array}\right)
V = M.row_space() W = M.right_kernel()
V;W
 Free module of degree 3 and rank 1 over Integer Ring Echelon basis matrix: [1 2 3] Free module of degree 3 and rank 2 over Integer Ring Echelon basis matrix: [ 1 1 -1] [ 0 3 -2] Free module of degree 3 and rank 1 over Integer Ring Echelon basis matrix: [1 2 3] Free module of degree 3 and rank 2 over Integer Ring Echelon basis matrix: [ 1 1 -1] [ 0 3 -2]

We just multiply every basis vector of the row space with every basis vector of the null space.

for v in V.basis(): for w in W.basis(): html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w)))
 $\left(1,\,2,\,3\right)\cdot \left(1,\,1,\,-1\right)=0$ $\left(1,\,2,\,3\right)\cdot \left(0,\,3,\,-2\right)=0$  
@interact def _(M = input_grid(4, 5, default=[[-1,-2,1,-1,0],[3,6,6,0,3],[3,6,4,1,3],[-2,-4,-5,0,-3]])): M = matrix(M) V = M.right_kernel() W = M.row_space() for v in V.basis(): for w in W.basis(): html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w)))

## Click to the left again to hide and once more to show the dynamic interactive window

Incidentally, the same is true of the transpose (the column space and kernel of the transpose):

@interact def _(M = input_grid(4, 5, default=[[-1,-2,1,-1,0],[3,6,6,0,3],[3,6,4,1,3],[-2,-4,-5,0,-3]])): M = matrix(M) V = M.left_kernel() W = M.column_space() html("Now we're doing it with the column space and <em>left</em> kernel") for v in V.basis(): for w in W.basis(): html("$%s\cdot %s=%s$"%(latex(v),latex(w),latex(v*w)))

## Click to the left again to hide and once more to show the dynamic interactive window

Warning!  There is a lot we are not saying about these complements - see Sage ticket #7522 for some discussion about this between those interested in undergraduate pedagogical situations and those warning of keeping things completely mathematically correct!

## Section 6.2

Let's do some orthogonal projection.

@interact def _(V = matrix([[1,3]]),W = matrix([[-2,3]])): V = vector(V.list()) W = vector(W.list()) ProjVW = (V*W)/(W*W)*W show(plot(ProjVW,color='green',zorder=5)+plot(V,color='green')+plot(W)+line([V.list(),ProjVW.list()],color='black',linestyle='--'),aspect_ratio=1)

## Section 6.3

Well, what about in more than two dimensions?

@interact def _(V = matrix([[1,3,2]]),W = matrix([[-2,3,2]])): V = vector(V.list()) W = vector(W.list()) ProjVW = (V*W)/(W*W)*W show(plot(ProjVW,color='green',zorder=5)+plot(V,color='green')+plot(W)+line([V.list(),ProjVW.list()],color='black',linestyle='--'),aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Note that the default setting here doesn't have the $W1$ and $W2$ being orthogonal!  You can't use the formula without doing that first.

var('u,v') @interact def _(V = matrix([[1,3,2]]),W1 = matrix([[-2,3,2]]),W2 = matrix([[-2,1,1]])): V = vector(V.list()) W1 = vector(W1.list()) W2 = vector(W2.list()) if W1*W2 != 0: print "you don't have orthogonal vectors" else: ProjVW = (V*W1)/(W1*W1)*W1+(V*W2)/(W2*W2)*W2 show(plot(ProjVW,color='green',zorder=5)+plot(V,color='blue',thickness=2)+parametric_plot3d(u*W1+v*W2, (u,-3,3), (v,-3,3),color='green',opacity=.7)+line([V.list(),ProjVW.list()],color='black',linestyle='--',thickness=5),aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

Of course, we can do all this for the polynomial or other inner products...

w1 = x w2 = 3/2*x^2-1/2 v = 1+x+x^2
F = integrate(w1*v,x,-1,1)/integrate(w1*w1,x,-1,1)*w1+integrate(w2*v,x,-1,1)/integrate(w2*w2,x,-1,1)*w2 F
 x^2 + x - 1/3 x^2 + x - 1/3

Here, $F$ is the projection and $v-F$ is the residual.

integrate((x^2+x-1/3)*(4/3),x,-1,1)
 0 0
integrate((v-F)*F,x,-1,1)
 0 0

Inner product is zero.  Yeah.  And the distance from $v$ to this subspace is...

sqrt(integrate((v-F)^2,x,-1,1))
 4/3*sqrt(2) 4/3*sqrt(2)

Huh.

## Section 6.4

So, how does one exactly get this orthogonal basis?  The answer is the Gram-Schmidt process!

#A = matrix([[1,1,0,1,0],[1,0,1,1,0]]) #A = matrix([[1,1,0,1,0],[1,0,1,1,0],[0,1,1,1,0]]) A = matrix([[1,1,0,1,0],[1,0,1,1,0],[0,1,1,1,0],[5,0,0,0,1]])
show(A.gram_schmidt()[0])
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} 1 & 1 & 0 & 1 & 0 \\ \frac{1}{3} & -\frac{2}{3} & 1 & \frac{1}{3} & 0 \\ -\frac{4}{5} & \frac{3}{5} & \frac{3}{5} & \frac{1}{5} & 0 \\ \frac{5}{7} & \frac{5}{7} & \frac{5}{7} & -\frac{10}{7} & 1 \end{array}\right)

Let's see what this looks like for three general vectors.

var('u,v') @interact def _(U1 = matrix([[1,3,2]]),U2 = matrix([[-2,3,2]]),U3 = matrix([[-2,1,1]])): U1 = vector(U1.list()) U2 = vector(U2.list()) U3 = vector(U3.list()) M = matrix([U1,U2,U3]) if M.det()==0: print "you don't have a basis!" else: M = M.gram_schmidt()[0] G = Graphics() G += plot(U1,thickness=2)+plot(U2,thickness=2)+plot(U3,thickness=2) G += parametric_plot3d(u*U1+v*U2, (u,-1,1), (v,-1,1),color='green',opacity=.6) V1,V2,V3 = M.rows() G += plot(V1,thickness=2,color='red')+plot(V2,thickness=2,color='red')+plot(V3,thickness=2,color='red') G += line([V2.list(),U2.list()],color='black',linestyle='--',thickness=2) G += line([V3.list(),U3.list()],color='black',linestyle='--',thickness=2) show(G,aspect_ratio=1)

## Click to the left again to hide and once more to show the dynamic interactive window

We can use a similar method to get a set of orthogonal polynomials - very useful in applications.

legendre_P?
 File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.py Type: Definition: legendre_P(n, x) Docstring: Returns the Legendre polynomial of the first kind for integers $$n > -1$$. REFERENCE: AS 22.5.35 page 779. EXAMPLES: sage: P. = QQ[] sage: legendre_P(2,t) 3/2*t^2 - 1/2 sage: legendre_P(3, 1.1) 1.67750000000000 sage: legendre_P(3, 1 + I) 7/2*I - 13/2 sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7])) [-179 242] [-484 547] sage: legendre_P(3, GF(11)(5)) 8  File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/functions/orthogonal_polys.py Type: Definition: legendre_P(n, x) Docstring: Returns the Legendre polynomial of the first kind for integers $$n > -1$$. REFERENCE: AS 22.5.35 page 779. EXAMPLES: sage: P. = QQ[] sage: legendre_P(2,t) 3/2*t^2 - 1/2 sage: legendre_P(3, 1.1) 1.67750000000000 sage: legendre_P(3, 1 + I) 7/2*I - 13/2 sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7])) [-179 242] [-484 547] sage: legendre_P(3, GF(11)(5)) 8 

## Section 6.5

Finally, let's fit some linear models.

var('a,b') model(x) = a+b*x Data = [(1,0),(2,1),(4,2),(5,3)] find_fit(Data,model)
 [a == -0.6000000000034891, b == 0.6999999999993456] [a == -0.6000000000034891, b == 0.6999999999993456]
M = find_fit(Data,model) points(Data)+plot(model.subs(a=M[0].rhs(),b=M[1].rhs()),(x,0,5),color='black')

Here's some real-life data about men's and women's record times in the 100-meter free in swimming.

var('a,b,c') model(gender,year) = a+b*year+c*gender Data = [(1,1905,65.8),(1,1908,65.6),(1,1910,62.8),(1,1912,61.6),(1,1918,61.4),(1,1920,60.4),(1,1922,58.6),(1,1924,57.4),(1,1934,56.8),(1,1935,56.6),(1,1936,56.4),(1,1944,55.9),(1,1947,55.8),(1,1948,55.4),(1,1955,54.8),(1,1957,54.6),(1,1961,53.6),(1,1964,52.9),(1,1967,52.6),(1,1968,52.2),(1,1970,51.9),(1,1972,51.22),(1,1975,50.59),(1,1976,49.44),(1,1981,49.36),(1,1985,49.24),(1,1986,48.74),(1,1988,48.42),(1,1994,48.21),(1,2000,48.18),(1,2000,47.84),(0,1908,95),(0,1910,86.6),(0,1911,84.6),(0,1912,78.8),(0,1915,76.2),(0,1920,73.6),(0,1923,72.8),(0,1924,72.2),(0,1926,70),(0,1929,69.4),(0,1930,68),(0,1931,66.6),(0,1933,66),(0,1934,65.4),(0,1936,64.6),(0,1956,62),(0,1958,61.2),(0,1960,60.2),(0,1962,59.5),(0,1964,58.9),(0,1972,58.5),(0,1973,57.54),(0,1974,56.96),(0,1976,55.65),(0,1978,55.41),(0,1980,54.79),(0,1986,54.73),(0,1992,54.48),(0,1994,54.01),(0,2000,53.77),(0,2004,53.52)]
M = find_fit(Data,model); M
 [a == 555.7167834308991, b == -0.2514636815413249, c == -9.797961450991004] [a == 555.7167834308991, b == -0.2514636815413249, c == -9.797961450991004]
model = model.subs(a=M[0].rhs(),b=M[1].rhs(),c=M[2].rhs()); model
 (gender, year) |--> -9.797961450991004*gender - 0.2514636815413249*year + 555.7167834308991 (gender, year) |--> -9.797961450991004*gender - 0.2514636815413249*year + 555.7167834308991
points([d[1:] for d in Data if d[0]==1])+points([d[1:] for d in Data if d[0]==0],color='red')+plot(model(gender=1),(year,1900,2010),color='blue',legend_label="modeling men's times" )+plot(model(gender=0),(year,1900,2010),color='red',legend_label="modeling women's times")