[Numpy-discussion] Total least squares problem
Nils Wagner
nwagner at mecha.uni-stuttgart.de
Wed Nov 14 04:44:03 CST 2001
Travis Oliphant schrieb:
>
> >
> > How do I solve a Total Least Squares problem in Numpy ?
> > A small example would be appreciated.
> >
> > The TLS problem assumes an overdetermined set of linear equations
> > AX = B, where both the data matrix A as well as the observation
> > matrix B are inaccurate:
>
> X, resids, rank, s = LinearAlgebra.linear_least_squares(A,B)
>
> -Travis
Travis,
There is a difference between classical least squares (Numpy)
and TLS (total least squares).
I am attaching a small example for illustration.
Nils
-------------- next part --------------
from Numeric import *
from LinearAlgebra import *
A = zeros((6,3),Float)
b = zeros((6,1),Float)
#
# Example by Van Huffel
# http://www.netlib.org/vanhuffel/dtls-doc
#
A[0,0] = 0.80010002
A[0,1] = 0.39985167
A[0,2] = 0.60005390
A[1,0] = 0.29996484
A[1,1] = 0.69990689
A[1,2] = 0.39997269
A[2,0] = 0.49994235
A[2,1] = 0.60003167
A[2,2] = 0.20012361
A[3,0] = 0.90013643
A[3,1] = 0.20016919
A[3,2] = 0.79995025
A[4,0] = 0.39998539
A[4,1] = 0.80006338
A[4,2] = 0.49985474
A[5,0] = 0.20002274
A[5,1] = 0.90007114
A[5,2] = 0.70009777
b[0] = 0.89999446
b[1] = 0.82997570
b[2] = 0.79011189
b[3] = 0.85002662
b[4] = 0.99016399
b[5] = 0.10299439
print 'Solution of an overdetermined system of linear equations A x = b'
print
print 'A'
print
print A
#
print 'b'
print
print b
#
x, resids, rank, s = linear_least_squares(A,b)
print
print 'Least squares solution (Numpy)'
print
print x
print
print 'Computed rank',rank
print
print 'Sum of the squared residuals', resids
print
print 'Singular values of A in descending order'
print
print s
#
xtls = zeros((3,1),Float)
#
# total least squares solution given by Van Huffel
# http://www.netlib.org/vanhuffel/dtls-doc
#
xtls[0] = 0.500254
xtls[1] = 0.800251
xtls[2] = 0.299492
print
print 'Total least squares solution'
print
print xtls
print
print 'Residuals of LS (Numpy)'
print
print matrixmultiply(A,x)-b
print
print 'Residuals of TLS'
print
print matrixmultiply(A,xtls)-b
print
#
# Least squares in Numpy A^\top A x = A^\top b
#
Atb = matrixmultiply(transpose(A),b)
AtA = matrixmultiply(transpose(A),A)
xls = solve_linear_equations(AtA,Atb)
print
print 'Least squares solution via normal equation'
print
print xls
More information about the Numpy-discussion
mailing list