Least Squares Approximation/Optimization - CSE SERVICES

Computational MethodsLeast SquaresApproximation/Optimization Manfred Huber 20111

Least Squaresn Least squares methods are aimed at findingapproximate solutions when no precise solutionexistsn n Find the solution that minimizes the residual error in thesystemLeast squares can be used to fit a model to noisydata points or to fit a simpler model to complexdatan Amounts to projecting higher dimensional data onto a lowerdimensional space. Manfred Huber 20112

Linear Least Squaresn Linear least squares attempts to find a leastsquares solution for an overdetermined linearsystem (i.e. a linear system described by an m x nmatrix A with! more equations than parameters).!Ax " bn Least squares minimizes the squared Eucliden norm of!!2!2the residualmin x! b " Ax min x! r 22!n For data fitting on m data points using a linearcombination of basis functions this corresponds tomnmin" % (y m # % " j j (x i )) 2i 1j 1! Manfred Huber 20113

Existence and Uniquenessn n Linear least squares problem always has a solutionSolution is unique if and only if A has full rank, i.e.rank(A) nn n If rank(A) n then A is rank-deficient and the solutionof the least squares problem is not uniqueIf solution is unique the residual vector can beexpressed through A2r 2 rT r (b " Ax)T (b " Ax) bT b " 2x T AT b x T AT Axn This is minimized if its derivative is 0 Manfred Huber 2011!"2AT b 2AT Ax 04

Normal Equationsn Optimization reduces to a n x n system of (linear)normal equationsTTA Ax A bn n Linear least squares can be found by solving this systemof linear equationsSolution can also be found through the Pseudo InverseAx " b # x A b!A (AT A) 1 ATn Condition number with respect to A can be expressed asin the case of solving linear equations! Manfred Huber 2011cond(A) A2A 25

Data Fittingn A common use for linear least squares solution is tofit a given type of function to noisy data pointsn nth order polynomial fit using monomial basis:#1 x1%! 1 x2A" %%# #% 1 x mn !# y1 &" x1n &% (n(" x2 ( ! % y2 (" %# ( # (% (n(" xm ' ym 'Solving the system of equations provides the best fit in terms ofthe Euclidian normTT Manfred Huber 2011A A" A y6

Condition Number andSensitivityn Sensitivity also depends on bn Influence of b can be expressed in terms of an anglebetween b and ycos(" ) n y2b2 Axb22Bound on the error in the solution due to perturbation inb can be expressed as!"xx Manfred Huber 201122"b 21# cond(A)cos( ) b 27

Condition Number andSensitivityn Bound on the error in the solution due toperturbation in A can be expressed as"xxn !n 222# ((cond(A)) tan( ) cond(A))"AA22For small residuals the condition number for least squaresis approximately cond(A).For large residuals the condition number can be square ofworse Manfred Huber 20118

Condition Number andSensitivityn Conditioning of normal equation solution isTcond(A A) (cond(A))n !n 2For large systems of equation the condition number ofthe formulation using normal equations (or thePseudoinverse) increases rapidlyMuch of the increased sensitivity is due to the need formultiplying A and AT in order to be able to apply asolution algorithm for the system of equationsn n Conditioning of the normal equations is potentially significantlyworse than the conditioning of the original systemAlgorithm is not very stable for large numbers of equations/datapoints Manfred Huber 20119

Augmented System Methodn Augmented system method transforms least squareproblem into an system of equation solvingproblem by adding equations and can be used toimprove the conditioningn Increase matrix size to a square (m n)x(m n) matrix byincluding the residual equationsr Ax bAT r 0n n "# I% T AA&# r & # b&(% ( % (0 ' x ' 0'Greater freedom to choose pivots and maintain stabilitySubstentially higher complexity for systems with m n due to theneed to solve a (n m)x(n m) system Manfred Huber 2011!10

QR Factorizationn n The Augmented System method addresses stabilitybut adds very high cost since it expands the matrixQR factorization changes the system matrix intosolvable form without computation of ATA andwithout expanding the matrix" R1,1 0 " A Q 0 0 "#0 Manfred Huber 2011R1,2 ! R1,n %'R2,2 ! R2,n '" # " '" R%'! 0 Rn,n ' Q '# 0&0 ! 0 ''" # " ''0 ! 0 &(QT Ax QT b11

QR Factorizationn n !The important part of Q for the solution consists of thefirst n rows since the others are multiplied by 0QT (Q1Q2 ) " Q1T Ax Rx Q1T bQR factorization factors the system matrix A into Q anR where R is upper triangularn As in LU Factorization, QT represents a sequence of solutionpreserving transformationsn n n In LU Factorization only identity has to be preserved which can be doneusing elimination matricesIn QR Factorization the least square characteristic has to be preservedwhich requires the use of orthogonal transformationsR is an upper triangular matrix that can be used for backwardsubstitution to compute the solution Manfred Huber 201112

Orthogonal Transformationsn Orthogonal transformations preserve the least squaresolutionQQT I2TTTTQv 2 (Qv) Qv v Q Qv v v vn !22For QR factorization we need to find a set of orthogonaltransformations that can transform A into an uppertriangular matrix R and that are numerically stablen Householder transformsn Givens rotationsn Gram-Schmidt orthogonalization Manfred Huber 201113

QR Factorization withHouseholder Transformationsn Householder transformations allow to zero all entriesbelow a chosen point in a vector an n Applying this consecutively for every column of the matrix A,choosing the diagonal element as the one below whicheverything is to be zeroed out we can construct RHouseholder transformation can be computed from a givenvector (vector of column values), a, setting all the values thatare not to be changed (entries above the diagonal) to 0vv TQ H I " 2 T , v a " #ei , # av vn !n 2The sign for β can be chosen to avoid cancellation (loss ofsignificance) in the computation of vH is orthogonal and symmetric: Manfred Huber 2011H H T H "114

QR Factorization withHouseholder Transformationsn In QR factorization with Householder transforms,successive columns are adjusted to upper triangular byzeroing all entries below the diagonal element.n Householder transform does not affect the entries in thecolumns to the left and the rows above the currently chosendiagonal element since it contains only 0s in the rows belowthe chosen diagonal term.n Applying the Householder transform only to the columns to the right savessignificant processing timeT!T!T &T#vava!vv !!vv !!!Ha j % I " 2 T ( a j a j " 2 T a j a j " 2v T j a j " 2 T j vv v'v vv vv v n Applying it to individual columns also eliminates the need to compute full!matrix H - only v a j " #e j is needed Manfred Huber 2011!15

QR Factorization Examplen Quadratic polynomial fit to 5 data pointsData : ("1,1),("0.5,0.5),(0,0),(0.5,0.5),(1,2)n !n !Design linear#1 x1%%1 x 2!A" %1 x 3%%1 x 4% 1 x 5system for data fitting# y1
x12 &(% (%y1x 22 (2(%%2( !x 3 " ) % y 3 ( * %1(% (%yx 42 (% 4(%1% (%(x 52 ' y5 ' 1# 1& 11 &(% ( 0.5 0.25(0.5(%!00 (" ) % 0 ((% (0.5 0.25(% 0.5((% (11 ' 2'Starting with the first column start eliminating entries belowthe diagonal using appropriate Householder transformschoosing the sign on β to avoid cancellation Manfred Huber 201116

QR Factorization Examplen Householder elimination by columnn 1st column:#1
& # 3.236&%(% ( %(101%(% ( %(!!v1 a1 " a1 2 e1 %1( 5 %0( % 1 (%(% ( %(101%% ((%% (( %%(( 1' 0' 1 'n !Negative sign on β because potentially cancelling term is 1#"2.236# "1.789&0"1.118&%(%(0"0.191"0.405"0.362%((! %H1 A % 00.309 "0.655( , H1 y % "0.862(%(%(00.809"0.405"0.362%%((%%((1.309 0.345 ' 0 1.138 ' Manfred Huber 201117

QR Factorization Examplen Householder elimination by columnn 2nd column:# 0 &# 0& # 0 &%(% ( %("0.1911"1.772%(% ( %(!!v 2 a2 " a2 2 e2 % 0.309 ( " 2.5 % 0( % 0.309 (%(% ( %(0.80900.809%%((%% (( %%(( 1.309 ' 0' 1.309 'n !Positive sign on β because possibly cancelling term is -0.191#"2.236#"1.789&0"1.118&%(%(01.58100.632%((! %H 2 H1 A % 00"0.725( , H 2 H1 y %"1.035(%(%(00"0.589"0.816%%((%%((00.047 ' 0 0.404 ' Manfred Huber 201118

QR Factorization Examplen Householder elimination by columnn 3rd column:# 0 &# 0& # 0 &%(% ( %(000%(% ( %(!!v 2 a3 " a3 2 e3 % "0.725( " 0.875 % 1( % "1.66 (%(% ( %("0.5890"0.589%%((%% (( %%(( 0.047 ' 0' 0.047 'n !Positive sign on β because possibly cancelling term is -0.725#"2.236#"1.789&0"1.118&%(%(01.58100.632%((! %H 3 H 2 H1 A % 000.935 ( , H 3 H 2 H1 y % 1.336 (%(%(0000.026%%((%%((00 ' 0 0.337 ' Manfred Huber 201119

QR Factorization Examplen Backward substitution with the upper triangular matrixyields the parameters and least squares fit secondorder polynomial#0.086&continuedExample,! %2(p(x) 0.086 0.4x 1.429x" % 0.4Resulting( curve and original data points are shown in graph%( 1.429 'Least Squares Data FittingExistence, Uniqueness, and ConditioningSolving Linear Least Squares ProblemsLeast SquaresData Fitting!! Manfred Huber 201120 interactive example

Other OrthogonalTransformationsn n Other transformations can be used for QR Factorizationn Givens rotationsn Gram-Schmidt OrthogonalizationHouseholder Transformations generally achieve thebest performance and stability tradeoffn Complexity of QR Factorization with Householdertransformations is approximately mn2-n3/3 multiplicationsn n Depending on the size of m (data points / equations) this is between thesame and two times the work of normal equationsConditioning of QR Factorization with Householdertransformations is optimal Manfred Huber 2011cond(A) r2(cond(A))221

QR Factorizationn Transformations for QR Factorization are numericallymore stable than elimination steps in LU Factorizationn Choice of sign in Householder transformations allows to avoidcancellation and thus instabilities in individual transformsn n Row Pivoting is usually not necessaryStability leads to QR factorization also frequently beingused instead of LU Factorization to solve nonsingularsystems of linear equationsn Increased complexity is traded off against stability (and thusprecision) of the solution Manfred Huber 201122

Nonlinear Least Squaresn As for equation solving, finding solutions for generalnonlinear systems requires iterative solutionsn Goal is to find the approximate solution with the smallestsquare residual, ρ, for the system of functions f(x)ri (x) bi " f i (x)1#(x) rT (x)r(x)2n At the minimum, the gradient of the square residual function,would be 0T"#(x) J r (x)r(x) 0!n Could use Multivariate Newton method to find the root of the derivativewhich would require second derivative, the Hessian of the square errorT Manfred Huber 2011!mH " (x) J r (x)J r (x) # ri (x)H ri (x)i 123

Gauss-Newton Methodn Multivariate Newton method would require solving thefollowing linear system in each iterationm#% T&J r (x)J r (x) " ri (x)H ri (x)(s )J rT (x)r(x)i 1 'n !n Requires frequent computation of the Hessian which isexpensive and reduces stabilityGauss-Newton avoids this by dropping second order termJ rT (x)J r (x)s "J T (x)r(x)n !This is best solved by converting this into the correspondingleast squares problem and using QR factorizationJ r (x)s " #r(x)n Once solved, Gauss-Newton operates like Multivariate Newton Manfred Huber 2011x t 1 x t s24

Gauss-Newton Methodn Gauss-Newton method replaces nonlinear least squaresproblem with a sequence of linear least squaresproblems that converge to solution of nonlinear systemn Converges if residual at the solution is not too largen n Large residual at solution can lead to large values in Hessian, potentiallyleading to slow convergence or, in extreme cases, non-convergenceIf it does not converge (large residual at solution) othermethods have to be usedn Levenberg-Marquardt method which uses an additional scalar parameter(and a separate, function-specific strategy to choose it) to modify step size(Jn (x)J r (x) µI) s "J (x)r(x)Tr# J r (x)' "r(x)'&)s * &)% µI ( % 0 (General optimization using the complete Hessian Manfred Huber 2011!Trn Significant increase in computational complexity25