213 - A14 - Least squares solutions

239 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