# 213 - A14 - Least squares solutions

## 815 days ago by Professor213

# This example is with three lines with no common point of intersection. # a = matrix(QQ,3,2,[1,-2,3,1,1,1]) # b = matrix(QQ,3,1,[-7,9,4]) # Here is an example which has two parallel lines a = matrix(QQ,3,2,[1,1,3,1,1,1]) b = matrix(QQ,3,1,[-7,9,4]) at = a.transpose() AtA = at*a print AtA soln = AtA.inverse()*at*b print soln print print soln.n() pt = (soln[0,0],soln[1,0])
 [11 5] [ 5 3] [ 21/4] [-27/4] [ 5.25000000000000] [-6.75000000000000] [11 5] [ 5 3] [ 21/4] [-27/4] [ 5.25000000000000] [-6.75000000000000]
var('x,y') # Let's compare the residuals using the Least Squares solution "pt" vs a host of other points. # and then graph it all. G = Graphics() x0 = pt[0] y0 = pt[1] others = (matrix([[x0+(random()-0.5)/2],[y0+(random()-0.5)/2]]) for k in range(10)) for x1 in others: print x1.transpose() print 'The residual using the above point is',norm(b-a*x1) print G += point((x1[0][0],x1[1][0]),color='green') print 'The residual using',(x0,y0),'is',norm(b-a*matrix([[x0],[y0]])) for k in range(3): print a[k,0],"*x+",a[k,1],"*y=",b[k,0] G += implicit_plot(a[k,0]*x+a[k,1]*y==b[k,0], (x,-10,10), (y,-10,10)) G += point(pt,color='red',size=20) G.show()
 [ 5.28974280667794 -6.96617375515508] The residual using the above point is 7.78277933371 [ 5.49426616838674 -6.82796164889984] The residual using the above point is 7.80923336152 [ 5.42321268151684 -6.80650849581778] The residual using the above point is 7.79369801268 [ 5.13815000572148 -6.92729273019398] The residual using the above point is 7.80578085002 [ 5.23292292782251 -6.54388592183457] The residual using the above point is 7.78430849081 [ 5.42698842192924 -6.70013584385187] The residual using the above point is 7.80642600693 [ 5.12780968181397 -6.97954759330883] The residual using the above point is 7.81682777346 [ 5.29266895524503 -6.75436721264346] The residual using the above point is 7.77934578296 [ 5.42138777350584 -6.69223699037310] The residual using the above point is 7.80590288755 [ 5.03066525593098 -6.65240601686821] The residual using the above point is 7.80023725748 The residual using (21/4, -27/4) is 7.77817459305 1 *x+ 1 *y= -7 3 *x+ 1 *y= 9 1 *x+ 1 *y= 4 [ 5.28974280667794 -6.96617375515508] The residual using the above point is 7.78277933371 [ 5.49426616838674 -6.82796164889984] The residual using the above point is 7.80923336152 [ 5.42321268151684 -6.80650849581778] The residual using the above point is 7.79369801268 [ 5.13815000572148 -6.92729273019398] The residual using the above point is 7.80578085002 [ 5.23292292782251 -6.54388592183457] The residual using the above point is 7.78430849081 [ 5.42698842192924 -6.70013584385187] The residual using the above point is 7.80642600693 [ 5.12780968181397 -6.97954759330883] The residual using the above point is 7.81682777346 [ 5.29266895524503 -6.75436721264346] The residual using the above point is 7.77934578296 [ 5.42138777350584 -6.69223699037310] The residual using the above point is 7.80590288755 [ 5.03066525593098 -6.65240601686821] The residual using the above point is 7.80023725748 The residual using (21/4, -27/4) is 7.77817459305 1 *x+ 1 *y= -7 3 *x+ 1 *y= 9 1 *x+ 1 *y= 4
# In 3d a = matrix(QQ,4,3,[1,3,5,1,1,0,1,1,2,1,3,3]) b = matrix(QQ,4,1,[3,5,7,-3]) at = a.transpose() AtA = at*a print AtA soln = AtA.inverse()*at*b print soln print print soln.n() pt = (soln[0,0],soln[1,0],soln[2,0])
 [ 4 8 10] [ 8 20 26] [10 26 38] [10] [-6] [ 2] [ 10.0000000000000] [-6.00000000000000] [ 2.00000000000000] [ 4 8 10] [ 8 20 26] [10 26 38] [10] [-6] [ 2] [ 10.0000000000000] [-6.00000000000000] [ 2.00000000000000]
var('x,y,z') k=0 G = Graphics() for k in range(4): print a[k,0],"*x+",a[k,1],"*y+",a[k,2],"*z=",b[k,0] G += implicit_plot3d(a[k,0]*x+a[k,1]*y+a[k,2]*z==b[k,0], (x,0,15), (y,-10,10), (z,-10,10),opacity=0.5) G += point3d(pt,color='red',size=20) G.show()
 1 *x+ 3 *y+ 5 *z= 3 1 *x+ 1 *y+ 0 *z= 5 1 *x+ 1 *y+ 2 *z= 7 1 *x+ 3 *y+ 3 *z= -3 1 *x+ 3 *y+ 5 *z= 3 1 *x+ 1 *y+ 0 *z= 5 1 *x+ 1 *y+ 2 *z= 7 1 *x+ 3 *y+ 3 *z= -3
# Let's let Sage compute the Least Squares solution to an underdetermined system # That is, a system which has infinitely many solutions...which is the "best" one. a = matrix(QQ,2,3,[1,3,5,1,1,0]) b = matrix(QQ,2,1,[3,5]) import numpy.linalg r = numpy.linalg.lstsq(a,b)[0] r
 array([[ 2.75925926], [ 2.24074074], [-1.2962963 ]]) array([[ 2.75925926], [ 2.24074074], [-1.2962963 ]])
x0=r[0][0] y0=r[1][0] z0=r[2][0] (x0,y0,z0)
 (2.75925925925926, 2.2407407407407427, -1.2962962962962981) (2.75925925925926, 2.2407407407407427, -1.2962962962962981)
print x0+3*y0+5*z0 print x0+y0
 3.0 5.0 3.0 5.0