Metoda najmniejszych kwadratów jako rozwiązanie normalne układu sprzecznego.



In [2]:

    
A=matrix(QQ,[[1,0],[0,1],[0,0]] ) 
b=vector( [1,3/4,2/3])
A.rank()









    Out[2]:





2



In [3]:

    
var('c1,c2')
show(A,'*',vector([c1,c2]).column(),"=",b.column())









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{array}\right) \verb|*| \left(\begin{array}{r} c_{1} \\ c_{2} \end{array}\right) \verb|=| \left(\begin{array}{r} 1 \\ \frac{3}{4} \\ \frac{2}{3} \end{array}\right)$



In [3]:

    
A\b









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-3681a29f48ae> in <module>()
----> 1 A * BackslashOperator() * b

/usr/lib/sagemath/local/lib/python2.7/site-packages/sage/misc/misc.pyc in __mul__(self, right)
   1107             (0.0, 0.5, 1.0, 1.5, 2.0)
   1108         """
-> 1109         return self.left._backslash_(right)
   1110 
   1111 

sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix._backslash_ (/usr/lib/sagemath//src/build/cythonized/sage/matrix/matrix2.c:4798)()

sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.solve_right (/usr/lib/sagemath//src/build/cythonized/sage/matrix/matrix2.c:6926)()

sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix._solve_right_general (/usr/lib/sagemath//src/build/cythonized/sage/matrix/matrix2.c:8296)()

ValueError: matrix equation has no solutions



In [27]:









    Out[27]:



In [13]:

    
A.columns()









    Out[13]:





[(1, 0, 0), (0, 1, 0)]



In [14]:

    
ars=sum([arrow3d((0,0,0),A.columns()[i] ) for i in range(2)])
barw=arrow3d((0,0,0),b,color='green')



In [15]:

    
plt=ars+barw
plt.show(aspect_ratio=[1,1,1], frame_aspect_ratio=[1,1,1],viewer='tachyon')



In [ ]:



In [16]:

    
lk=A.left_kernel().basis()[0]
lk









    Out[16]:





(0, 0, 1)



In [ ]:



In [18]:

    
lka=arrow3d((0,0,0),lk,color='red')



In [19]:

    
plt=lka+ars+barw
plt.show(aspect_ratio=[1,1,1], frame_aspect_ratio=[1,1,1],viewer='tachyon')



In [ ]:



In [ ]:

Zamiast:

$$ A x = b$$

rozwiązujemy

$$ A^T A x = A^T b.$$



In [11]:

    
sol_lsq=A.transpose()*A\A.transpose()*b



In [12]:

    
A*sol_lsq









    Out[12]:





(1, 3/4, 0)



In [13]:

    
sol_lsq_arw=arrow3d((0,0,0),A*sol_lsq,color='yellow')



In [14]:

    
plt=lka+ars+barw+sol_lsq_arw
plt.show(aspect_ratio=[1,1,1], frame_aspect_ratio=[1,1,1],viewer='tachyon')

Inaczej mówiąc, niech błąd:

$$ r=b-A x$$

leży w lewym jądrze A:

$$A^T r =A^T( b-Ax) = A^T b-A^TAx = 0.$$

$$ A^T A x = A^T b.$$



In [33]:









    Out[33]:

sagecell 3d



In [21]:

    
#A=matrix(QQ,[[1,3],[2,1],[1,2]] ) 
#A=matrix(QQ,[[1,0],[0,1],[0,0]] ) 
A=random_matrix(QQ,3,2)
b=vector( [1,2/3,3/4])
print A.rank()


ars=sum([arrow3d((0,0,0),A.columns()[i] ) for i in range(2)])
barw=arrow3d((0,0,0),b,color='green')

lk=A.left_kernel().basis()[0]
lka=arrow3d((0,0,0),lk,color='red')+arrow3d((0,0,0),-lk,color='red')

sol_lsq=A* (A.transpose()*A\A.transpose()*b)
sol_lsq_arw=arrow3d((0,0,0),sol_lsq,color='yellow',viewer='tachyon')

l3d=line3d([ b,sol_lsq],thickness=3,color='gray')

plt=lka+ars+barw+sol_lsq_arw+l3d
plt.show(aspect_ratio=[1,1,1], frame_aspect_ratio=[1,1,1],viewer='tachyon')
show(A,b)









    



2






    












    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & -1 \\ 0 & 1 \\ 0 & 0 \end{array}\right) \left(1,\,\frac{2}{3},\,\frac{3}{4}\right)$



In [ ]:



In [31]:









    Out[31]:



In [32]:









    Out[32]:



In [34]:









    Out[34]: