[Numpy-discussion] Total least squares problem

[Nils Wagner]

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