> # HW#5 Key > > with(linalg): > grade:=0; Warning, the protected names norm and trace have been redefined and unprotected grade := 0 > > #1. Let the following be a list of points in R^3 > > points := [[0,0,6.98], > [1,-1,18.01], > [1,1,10.01], > [-1,1,3.97], > [-1,-1,8.04], > [1,2,6.02], > [-2,3,7.98], > [2,3,12.03], > [2,4,7.04], > [2,2,16.99]]; points := [[0, 0, 6.98], [1, -1, 18.01], [1, 1, 10.01], [-1, 1, 3.97], [-1, -1, 8.04], [1, 2, 6.02], [-2, 3, 7.98], [2, 3, 12.03], [2, 4, 7.04], [2, 2, 16.99]] > # We'll assume that these points all lie on the plane z=a+b*x+c*y. > We'll want to find a, b, and c using Least Squares. > #a. Find a matrix A and a vector b such that if x=[a,b,c]^T then a > solutions to A.x=b would yield the coefficients of the plane. > > #For every point, the following should be true: z=a+b*x+c*y or > z=[1,x,y].[a,b,c]. So the rows of the matrix A should have the form > [1,x,y] > n:=nops(points): > A1:=matrix(n,3,0): > b1:=vector(n,0): > for i from 1 to n do > A1[i,1] := 1; A1[i,2]:=points[i][1]; A1[i,3]:=points[i][2]; > b1[i]:=points[i][3]; > end do: > evalm(A1); evalm(b1); [1 0 0] [ ] [1 1 -1] [ ] [1 1 1] [ ] [1 -1 1] [ ] [1 -1 -1] [ ] [1 1 2] [ ] [1 -2 3] [ ] [1 2 3] [ ] [1 2 4] [ ] [1 2 2] [6.98, 18.01, 10.01, 3.97, 8.04, 6.02, 7.98, 12.03, 7.04, 16.99] > grade:=grade+3; grade := 3 > # b. Do a least squares fit to find x^* and b^* such that A.x^* = b^*. > > x_star:=leastsqrs(A1,b1); x_star := [10.17889850, 2.000905780, -1.051679563] > b_star:=evalm(A1&*x_star); b_star := [10.17889850, 13.23148384, 11.12812472, 7.126313157, 9.229672283, 10.07644515, 3.022048251, 11.02567137, 9.973991808, 12.07735093] > grade:=grade+3; grade := 6 > > # c. How close is the least squares fit? > norm(b1-b_star,2); 10.97643575 > # Not such a hot fit. > grade:=grade+2; grade := 10 > > #2. What if the points from problem 1 lie on a different type of > algebraic surface? Suppose they lie on a quadratic surface: > z=a+b*x+c*y+d*xy+e*x^2. Repeat problem 1 for this matrix, though this > time x=[a,b,c,d,e]^T. Which of these is a better fit for the data? > Why? > > #a. > n:=nops(points): > A2:=matrix(n,5,0): > b2:=b1; # Same z points. > for i from 1 to n do > A2[i,1] := 1; A2[i,2]:=points[i][1]; A2[i,3]:=points[i][2]; > A2[i,4]:=A2[i,2]*A2[i,3]; A2[i,5]:=A2[i,2]^2; > end do: > evalm(A2); evalm(b2); b2 := b1 [1 0 0 0 0] [ ] [1 1 -1 -1 1] [ ] [1 1 1 1 1] [ ] [1 -1 1 -1 1] [ ] [1 -1 -1 1 1] [ ] [1 1 2 2 1] [ ] [1 -2 3 -6 4] [ ] [1 2 3 6 4] [ ] [1 2 4 8 4] [ ] [1 2 2 4 4] [6.98, 18.01, 10.01, 3.97, 8.04, 6.02, 7.98, 12.03, 7.04, 16.99] > grade:=grade+3; grade := 13 > #b. > x2_star:=leastsqrs(A2,b2); x2_star := [7.001836981, 3.994503624, -3.005128185, -.9936412651, 3.002946874] > b2_star:=evalm(A2&*x2_star); b2_star := [7.001836981, 17.99805692, 10.00051802, 3.998793311, 8.021767151, 6.001748574, 7.971080271, 12.02539958, 7.032988869, 17.01781030] > > #c. > norm(b2-b_star,2); > norm(b2-b2_star,2); 10.97643575 .05592225951 > # Compairing the norms, it seems that these points are better fit to a > quadratic surface. > grade:=grade+3; grade := 16 > > # 3. Given x, and y, in R^3. Determine if the following are inner > products in R^3. Explain why or why not. > > x:=vector(3): y:=vector(3): > > # Some sample vectors: > u := vector([1,1,0]): > v := vector([-1,0,1]): > e1 := vector([1,0,0]): e2 := vector([0,1,0]): e3 := vector([0,0,1]): > > > #a. = x[1]^2*y[1] + x[2]*y[2] +x[3]*y[3]; > ipa := (x,y) -> x[1]^2*y[1] + x[2]*y[2] +x[3]*y[3]; 2 ipa := (x, y) -> x[1] y[1] + x[2] y[2] + x[3] y[3] > ipa(x,y); 2 x[1] y[1] + x[2] y[2] + x[3] y[3] > ipa(evalm(u+v),v) - (ipa(u,v) + ipa(v,v)); # Check bilinearity > ipa(evalm(2*u),v) - 2*ipa(u,v); 2 -2 > # This is not bilinear. Thus not an innerproduct > grade:=grade+2; grade := 18 > #b. = x[1]^2+y[2]^2+x[3]*y[3] > ipb := (x,y) -> x[1]^2+y[2]^2+x[3]*y[3]; 2 2 ipb := (x, y) -> x[1] + y[2] + x[3] y[3] > ipb(e1,e2) - ipb(e2,e1); # Check symmetry 2 > # Not symmetric, so not an innerproduct. > grade:=grade+2; grade := 20 > > #c. =|x1|+|y1|+|x2|+|y2|+|x3|+|y3| > ipc:=(x,y)->abs(x[1])+abs(x[2])+abs(x[3])+abs(y[1])+abs(y[2])+abs(y[3] > ); ipc := (x, y) -> | x[1] | + | x[2] | + | x[3] | + | y[1] | + | y[2] | + | y[3] | > ipc(evalm(-1*u),v) - (-1*ipc(u,v)); # Check bilinearity 8 > # Not bilinear, so not an innerproduct > grade:=grade+2; grade := 22 > > #d. = 2*x1*y1 - .00001*x2*y2 + 1000000*x3*y3 > ipd:=(x,y)->2*x[1]*y[1] - .00001*x[2]*y[2] + 1000000*x[3]*y[3]; ipd := (x, y) -> 2 x[1] y[1] - .00001 x[2] y[2] + 1000000 x[3] y[3] > ipd(e2,e2); # Check if this is positive -.00001 > # can be negative, so this is not an inner product. > grade:=grade+2; grade := 24 > #e. = 2*x1*y1 + .00001*x2*y2 + 1000000*x3*y3 > ipe:=(x,y) -> 2*x[1]*y[1] + .00001*x[2]*y[2] + 1000000*x[3]*y[3]; > z:=vector(3): ipe := (x, y) -> 2 x[1] y[1] + .00001 x[2] y[2] + 1000000 x[3] y[3] > ipe(x,y)-ipe(y,x); # It's symmetric 0. > expand(ipe(evalm(x+y),z))-(ipe(x,z)+ipe(y,z)); # It's bilinear > expand(ipe(evalm(a*x),z))-expand(a*ipe(x,z)); 0. 0. > ipe(x,x); # This is always positive. 2 2 2 2 x[1] + .00001 x[2] + 1000000 x[3] > grade:=grade+2; grade := 24 > > #4. Let u, v, w, z, and b be as follows: > u4:=vector([1,0,-1,0,3]); > v4:=vector([-1,0,2,-1,0]); > w4:=vector([4,-3,1,0,0]); > z4:=vector([6,-3,-2,1,3]); > b4:=vector([5,7,-4,3,-2]); u4 := [1, 0, -1, 0, 3] v4 := [-1, 0, 2, -1, 0] w4 := [4, -3, 1, 0, 0] z4 := [6, -3, -2, 1, 3] b4 := [5, 7, -4, 3, -2] > #a. Using the Gram-Schmidt process, compute an orthogonal basis for > the span(u,v,w,z) > new_basis:=GramSchmidt([u4,v4,w4,z4]); -8 19 new_basis := [[1, 0, -1, 0, 3], [--, 0, --, -1, 9/11], 11 11 203 -13 -12 [---, -3, 5/3, ---, ---]] 57 57 19 > > grade:=grade+3; grade := 27 > #b. Using the basis computed in a. find the projection of b onto the > span. > > #Using the new basis, which is orthogonal, this is easy as pie. > > b_hat := dotprod(new_basis[1],b4)/dotprod(new_basis[1],new_basis[1]) * > new_basis[1] > + dotprod(new_basis[2],b4)/dotprod(new_basis[2],new_basis[2]) * > new_basis[2] > + dotprod(new_basis[3],b4)/dotprod(new_basis[3],new_basis[3]) * > new_basis[3]; 1529 1587 -1691 4281 -477 b_hat := [----, ----, -----, ----, ----] 1420 1420 284 1420 355 > # If you did problem #3a in Homework #4 correctly, this is the same > thing you computed there. > grade:=grade+3; grade := 30 > > > > >