Hi, 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: Nils Reference: R.D.Fierro, G.H. Golub, P.C. Hansen, D.P.O'Leary, Regularization by truncated total least squares, SIAM J. Sci. Comput. Vol.18(4) 1997 pp. 1223-1241
Travis Oliphant schrieb:
Travis, There is a difference between classical least squares (Numpy) and TLS (total least squares). I am attaching a small example for illustration. Nils 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
Nils Wagner <nwagner@mecha.uni-stuttgart.de> writes:
There is a difference between classical least squares (Numpy) and TLS (total least squares).
Algorithmically speaking it is even a very different problem. I'd say the only reasonable (i.e. efficient) solution for NumPy is to implement the TLS algorithm in a C subroutine calling LAPACK routines for SVD etc. Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hinsen@cnrs-orleans.fr Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais -------------------------------------------------------------------------------
Konrad Hinsen schrieb:
There are two Fortran implementations of the TLS algorithm already available via http://www.netlib.org/vanhuffel/ . Moreover there is a tool called f2py that generates Python C/API modules for wrapping Fortran 77/90/95 codes to Python. Unfortunately I am not very familar with this tool. Therefore I need some advice for this. Thanks in advance Nils
Travis Oliphant schrieb:
Travis, There is a difference between classical least squares (Numpy) and TLS (total least squares). I am attaching a small example for illustration. Nils 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
Nils Wagner <nwagner@mecha.uni-stuttgart.de> writes:
There is a difference between classical least squares (Numpy) and TLS (total least squares).
Algorithmically speaking it is even a very different problem. I'd say the only reasonable (i.e. efficient) solution for NumPy is to implement the TLS algorithm in a C subroutine calling LAPACK routines for SVD etc. Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hinsen@cnrs-orleans.fr Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais -------------------------------------------------------------------------------
Konrad Hinsen schrieb:
There are two Fortran implementations of the TLS algorithm already available via http://www.netlib.org/vanhuffel/ . Moreover there is a tool called f2py that generates Python C/API modules for wrapping Fortran 77/90/95 codes to Python. Unfortunately I am not very familar with this tool. Therefore I need some advice for this. Thanks in advance Nils
participants (3)
-
Konrad Hinsen -
Nils Wagner -
Travis Oliphant