# 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])
show(a)
b = matrix(QQ,3,1,[-7,9,4])
show(b)
at = a.transpose()
AtA = at*a
show(AtA)
soln = AtA.inverse()*at*b
show(soln)
print
r = a*soln-b
show(r)
print soln.n()
pt = (soln[0,0],soln[1,0])
show(r.norm())
[ 5.25000000000000] [-6.75000000000000] [ 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 |
|
|