Näherungslösungen von überbestimmten Gleichungssystemen

Gegeben: $A \cdot x = b$ wobei $A \in M(m \times n), m > n$
Gesucht: Näherungslösung $\bar{x}$ sodass $||b - A\bar{x}|| \leq || b - Ax||$



In [1]:

    
# input: equation system Ax = b
A = matrix([
        [1, 2],
        [2, 3],
        [3, 4]])
b = vector([1, 
            2,
            0])
#

Mit "einfacher" Pseudo-Inverse

$$\Sigma^{\#} = (A^T \cdot A)^{-1} \cdot A^T$$



In [2]:

    
Apinv = (A.T * A)^-1 * A.T
show(LatexExpr(r"A{T} \cdot A = "), A.T * A)
show(LatexExpr(r"(A{T} \cdot A)^{-1} = "), (A.T * A)^-1)
show(LatexExpr(r"A^{\#} = (A{T} \cdot A)^{-1} \cdot A^{T} = "), Apinv)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}A{T} \cdot A = \left(\begin{array}{rr} 14 & 20 \\ 20 & 29 \end{array}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}(A{T} \cdot A)^{-1} = \left(\begin{array}{rr} \frac{29}{6} & -\frac{10}{3} \\ -\frac{10}{3} & \frac{7}{3} \end{array}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}A^{\#} = (A{T} \cdot A)^{-1} \cdot A^{T} = \left(\begin{array}{rrr} -\frac{11}{6} & -\frac{1}{3} & \frac{7}{6} \\ \frac{4}{3} & \frac{1}{3} & -\frac{2}{3} \end{array}\right)$



In [3]:

    
x = Apinv * b
show(LatexExpr(r"\bar{x} ="), x)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\bar{x} = \left(-\frac{5}{2},\,2\right)$



In [4]:

    
show(LatexExpr(r"A \cdot \bar{x} ="), A*x)
show(LatexExpr(r"b = "), b)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}A \cdot \bar{x} = \left(\frac{3}{2},\,1,\,\frac{1}{2}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}b = \left(1,\,2,\,0\right)$



In [ ]:

Mit Pseudoinverse ermittelt mit SVD

Singulärwertzerlegung $A = U \cdot \Sigma V^{T}$

Die Pseudoinverse ist definiert als: $$ A^{\#} = V \cdot \Sigma^{\#} \cdot U^{T} $$ mit $\Sigma^{\#} = \begin{pmatrix}D^{-1} & 0 \\ 0 & 0 \\ \end{pmatrix} \in M(n \times m)$



In [5]:

    
U, S, V  = matrix(RDF, A).SVD()
show(U,S,V)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -0.338098165838 & 0.847952217752 & 0.408248290464 \\ -0.550649318286 & 0.173547289246 & -0.816496580928 \\ -0.763200470734 & -0.50085763926 & 0.408248290464 \end{array}\right) \left(\begin{array}{rr} 6.54675563644 & 0.0 \\ 0.0 & 0.37415322624 \\ 0.0 & 0.0 \end{array}\right) \left(\begin{array}{rr} -0.569594837763 & -0.821925617556 \\ -0.821925617556 & 0.569594837763 \end{array}\right)$



In [6]:

    
Si = matrix(RDF, A.ncols(), A.nrows())
for i in range(S.ncols()):
    Si[i,i] = S[i,i]^-1 if S[i,i] != 0 else 0.0
Apinv = V *  Si * U.T
show(Apinv)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -1.83333333333 & -0.333333333333 & 1.16666666667 \\ 1.33333333333 & 0.333333333333 & -0.666666666667 \end{array}\right)$



In [7]:

    
x = Apinv * b
show(LatexExpr(r"\bar{x} ="), x)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\bar{x} = \left(-2.5,\,2.0\right)$



In [8]:

    
show(LatexExpr(r"A \cdot \bar{x} ="),A*x)
show(LatexExpr(r"b = "), b)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}A \cdot \bar{x} = \left(1.5,\,1.0,\,0.5\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}b = \left(1,\,2,\,0\right)$



In [ ]:

Mit QR Zerlegung

QR Zerlegung $A = QR$

Näherungslösung wird berechnet durch lösen von: $$R \bar{x} = Q^T b$$



In [9]:

    
try:
    # try rational field first
    Q,R = A.QR()
except TypeError:
    # if it doesn't work, use real double field
    Q,R = matrix(RDF, A).QR()
    
show(LatexExpr(r"Q ="), Q)
show(LatexExpr(r"R ="), R)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}Q = \left(\begin{array}{rrr} -0.267261241912 & 0.872871560944 & 0.408248290464 \\ -0.534522483825 & 0.218217890236 & -0.816496580928 \\ -0.801783725737 & -0.436435780472 & 0.408248290464 \end{array}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}R = \left(\begin{array}{rr} -3.74165738677 & -5.34522483825 \\ 0.0 & 0.654653670708 \\ 0.0 & 0.0 \end{array}\right)$



In [10]:

    
Xsym = vector([var("x{}".format(i)) for i in range(A.ncols())])
RX = R * Xsym
QTb = Q.T * b
show(LatexExpr(r"R x = "), RX)
show(LatexExpr(r"Q^T b = "),QTb)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}R x = \left(-3.74165738677 \, x_{0} - 5.34522483825 \, x_{1},\,0.654653670708 \, x_{1},\,0.0\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}Q^T b = \left(-1.33630620956,\,1.30930734142,\,-1.22474487139\right)$



In [11]:

    
sol = solve([RX[i] == QTb[i] for i in range(A.nrows())], *Xsym)
if sol:
    for s in sol:
        show(s)
else:
    print "No solution found"









    



No solution found



In [ ]: