# The NonlinearProgramming Package # by Jason Schattman, Waterloo Maple, Inc., October 16, 2000 # jschattman@maplesoft.com # restart; # The Package Code NonlinearProgramming := module() export Optimize, PrimalDualLogBarrier, UnconstrainedNewton; local QuadraticIsConvex, FunctionIsQuadratic, FunctionIsLinear, SetIsLinear, IsPositiveVector, zerorhs, EvaluateVector, EvaluateDerivative, EvaluateDerivativeSum, VectorSubtractFast, MatrixSubtractFast, VectorAddFast, MatrixAddFast, ScaledVectorAdd, ScaledMatrixAdd, MultiplyMM, MultiplyMV, MultiplyMS, printConstraints, MultiplyVS, LinearSolveFast, TransposeFast, MatrixColumn, VectorNormFast, DiagMatrix, SubVectorFast, DotProd, IsPosDefinite, DimensionFast, ZeroVectorFast, ScaledVectorAddInPlace, ScaledMatrixAddInPlace, AddBoundConstraints, GetConstraintFunctions, ConvertConstraints, PickAlgorithm, OptimizeSimplex, OptimizeExplicitFONOC, OptimizeInteriorPointMethod, ShouldStop, InitialBarrierParameter, ComputeBarrierParameter, ComputeDirectionXW, ComputeDirectionY, ComputeDirectionW, FeasibleStepSize, MeritStepSize, ConstrainedArmijoStepSize, ComplexArmijoStepSize, SetMeritPenaltyParameter, ShouldTerminate; option package; QuadraticIsConvex := proc( hess ) option inline; LinearAlgebra:-LA_Main:-IsDefinite( hess, 'query=positive_def' ); end proc: IsPosDefinite := proc( M ) option inline; LinearAlgebra:-LA_Main:-IsDefinite( M, 'query=positive_def' ); end proc: FunctionIsQuadratic := proc( f, n ) option inline; type( f, quadratic( [seq( x[i], i=1..n )] ) ); end proc: FunctionIsLinear := proc( f, n ) option inline; type( f, linear( [seq( x[i], i=1..n )] ) ); end proc: SetIsLinear := proc( functionSet, n ) local i; for i from 1 to nops( functionSet ) do if not FunctionIsLinear( functionSet[i], n ) then return false; end if; end do; true; end proc: IsPositiveVector := proc( v, dim ) local i; for i from 1 to dim do if v[i] < 0 then return false; end if; end do; true; end proc: zerorhs := proc( inequal ) # converts the inequality to form h( x )>=0 and returns h( x ) option inline; -`-`( op( inequal ) ); end proc: EvaluateVector := proc( vector_expr, arguments, dim ) option inline; subs( {seq( x[i]=arguments[i], i=1..dim )}, vector_expr ); end proc: EvaluateDerivative := proc( derivative_expr, arguments, dim ) local expression; if type( derivative_expr, 'Vector' ) then expression := TransposeFast( derivative_expr ); else expression := derivative_expr; end if; subs( {seq( x[i]=arguments[i], i=1..dim )}, expression ); end proc: EvaluateDerivativeSum := proc( derivative_expr, arguments1, arguments2, dim1, dim2 ) local expression; if type( derivative_expr, 'Vector' ) then expression := TransposeFast( derivative_expr ); else expression := derivative_expr; end if; subs( {seq( x[i]=arguments1[i], i=1..dim1 ), seq( k[i]=arguments2[i], i=1..dim2 )}, expression ); end proc: VectorSubtractFast := proc( v1, v2 ) option inline; LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, 1, -1, 'inplace=false', 'outputoptions'=[] ); end proc: MatrixSubtractFast := proc( M1, M2 ) option inline; LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, 1, -1, 'inplace=false', 'outputoptions'=[] ); end proc: VectorAddFast := proc( v1, v2 ) option inline; LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, 1, 1, 'inplace=false', 'outputoptions'=[] ); end proc: MatrixAddFast := proc( M1, M2 ) option inline; LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, 1, 1, 'inplace=false', 'outputoptions'=[] ); end proc: ScaledVectorAdd := proc( v1, v2, scale1, scale2 ) option inline; LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, scale1, scale2, 'inplace=false', 'outputoptions'=[] ); end proc: ScaledMatrixAdd := proc( M1, M2, scale1, scale2 ) option inline; LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, scale1, scale2, 'inplace=false', 'outputoptions'=[] ); end proc: MultiplyMM := proc( M1, M2 ) option inline; LinearAlgebra:-LA_Main:-MatrixMatrixMultiply( M1, M2, 'inplace=false', 'outputoptions'=[] ); end proc: MultiplyMV := proc( M1, M2 ) option inline; LinearAlgebra:-LA_Main:-MatrixVectorMultiply( M1, M2, 'outputoptions'=[] ); end proc: MultiplyMS := proc( M1, s ) option inline; LinearAlgebra:-LA_Main:-MatrixScalarMultiply( M1, s, 'inplace=false', 'outputoptions'=[] ); end proc: MultiplyVS := proc( v1, s ) option inline; LinearAlgebra:-LA_Main:-VectorScalarMultiply( v1, s, 'inplace=false', 'outputoptions'=[] ); end proc: LinearSolveFast := proc( M, v ) option inline; LinearAlgebra:-LA_Main:-LinearSolve(M, v, 'method=Cholesky', 'free=s', 'conjugate=false', 'inplace=false', 'outputoptions'=[], 'methodoptions'=[] ); end proc: TransposeFast := proc( M ) option inline; LinearAlgebra:-LA_Main:-Transpose( M, 'inplace=false', 'outputoptions'=[] ); end proc: MatrixColumn := proc( M, c ) option inline; LinearAlgebra:-LA_Main:-Column( M, c, 'outputoptions'=[] ); end proc: VectorNormFast := proc( v ) option inline; LinearAlgebra:-LA_Main:-VectorNorm( v, 2 ); end proc: DiagMatrix := proc( veclist, n) option inline; LinearAlgebra:-LA_Main:-DiagonalMatrix( veclist, n, n,'outputoptions'=[datatype=float[8]] ); end proc: SubVectorFast := proc( v, index ) option inline; LinearAlgebra:-LA_Main:-SubVector( v, index, 'outputoptions'=[] ); end proc: DotProd := proc( v1, v2 ) option inline; LinearAlgebra:-LA_Main:-DotProduct( v1, v2, 'conjugate=false' ); end proc: DimensionFast := proc( M ) option inline; LinearAlgebra:-LA_Main:-Dimension( M ); end proc: ZeroVectorFast := proc( dim ) option inline; LinearAlgebra:-LA_Main:-ZeroVector[column]( dim, 'compact=false', 'outputoptions'=[] ); end proc: ScaledVectorAddInPlace := proc( v1, v2, scale1, scale2 ) option inline; LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, scale1, scale2, 'inplace=true', 'outputoptions'=[] ); end proc: ScaledMatrixAddInPlace := proc( M1, M2, scale1, scale2 ) option inline; LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, scale1, scale2, 'inplace=true', 'outputoptions'=[] ); end proc: AddBoundConstraints := proc( constraintSet::set, variableBounding, numVars::integer ) if variableBounding = 'free' then constraintSet; elif variableBounding = 'nonnegative' then constraintSet union {seq(x[i]>=0, i = 1..numVars)}; else error "Variable bounds must have value 'free' or 'nonnegative' but received %2", variableBounding; end if; end proc: ConvertConstraints := proc( constraintSet::set, variableBounding, numVars::integer ) local constraintsWithBounds, constraintFunctions; constraintsWithBounds := AddBoundConstraints( constraintSet, variableBounding, numVars ); constraintFunctions := GetConstraintFunctions( constraintsWithBounds ); Vector( convert( constraintFunctions, list ) ); end proc: GetConstraintFunctions := proc( constraintSet::set ) map( zerorhs, constraintSet ); end proc: printConstraints := proc ( constraints, m ) [ seq(constraints[i]>=0, i=1..m) ]; end proc: Optimize := proc(function, objective, num_variables, constraint_set, variable_bounds, x_start) local obj_function, constraint_functions, num_constraints, all_constraints; if objective='max' then obj_function := -1 * function; userinfo( 2, 'Optimize', 'Maximize' ); elif objective='min' then obj_function := function; userinfo( 2, 'Optimize', 'Minimize' ); else error "expecting the %-1 argument to be either `%2` or `%3`", 2, 'min', 'max'; end if; userinfo( 1, 'Optimize', function ); PickAlgorithm( obj_function, num_variables, constraint_set, variable_bounds, x_start ); end proc: PickAlgorithm := proc( obj_function, numVars, constraintSet, variableBounding, x_start ) local n, all_constraints, constraintFunctions, constraintStructure, hess; n := numVars; if nops( constraintSet ) = 0 and variableBounding = 'free' then constraintStructure := unconstrained; else constraintStructure := constrained; end if; if SetIsLinear( GetConstraintFunctions( constraintSet ), n ) then if FunctionIsLinear( obj_function, n ) then OptimizeSimplex( obj_function, constraintSet, variableBounding, n ); elif FunctionIsQuadratic( obj_function, n ) then hess := Matrix( hessian( obj_function, [seq( x[i], i=1..n )] ) ); if QuadraticIsConvex( hess ) then if constraintStructure = unconstrained then OptimizeExplicitFONOC( obj_function, hess, n ); else OptimizeInteriorPointMethod( obj_function, n, constraintSet, variableBounding, x_start, convex, constrained ); end if; else OptimizeInteriorPointMethod( obj_function, n, constraintSet, variableBounding, x_start, nonconvex, constraintStructure ); end if; else OptimizeInteriorPointMethod( obj_function, n, constraintSet, variableBounding, x_start, nonconvex, constraintStructure ); end if; else OptimizeInteriorPointMethod( obj_function, n, constraintSet, variableBounding, x_start, nonconvex, constraintStructure ); end if; end proc: OptimizeSimplex := proc( objFunction, constraintSet, variableBounding, numVars ) local constraintsWithBounds; userinfo( 2, 'Optimize', `Linear problem: solving with simplex method` ); constraintsWithBounds := AddBoundConstraints( constraintSet, variableBounding, numVars ); simplex['minimize']( objFunction, constraintsWithBounds ); end proc: OptimizeExplicitFONOC := proc( obj_function, hess, n ) local gradient, gradient0, soln; userinfo( 2, 'Optimize', `Convex quadratic problem: solving explicitly with first-order conditions` ); gradient := Vector( grad( obj_function, [seq( x[i], i=1..n )] ) ); gradient0 := EvaluateDerivative( gradient, ZeroVectorFast( n ), n ); userinfo( 2, 'OptimizeExplicitFONOC', `Global optimum found at` ); soln := LinearSolveFast( hess, MultiplyVS( gradient0, -1 ) ); { seq( x[i] = soln[i], i = 1..n) }; end proc: OptimizeInteriorPointMethod := proc( obj_function, numVars, constraintSet, variableBounding, x_start, convexityType, constraintStructure ) if constraintStructure = constrained then userinfo( 2, 'Optimize', `Constrained`, convexityType, `problem: solving with`, convexityType, `primal-dual log-barrier algorithm`); PrimalDualLogBarrier( obj_function, numVars, constraintSet, variableBounding, x_start, convexityType); else userinfo( 2, 'Optimize', `Unconstrained`, convexityType, `problem: solving with`, convexityType, `Newton's Method`); UnconstrainedNewton( obj_function, numVars, x_start, convexityType, float[8] ); end if; end proc: PrimalDualLogBarrier := proc(obj_function, numVars, constraintSet, variableBounding, x_start, convexityType) local x_curr, w_slack, y_mult, n, m, fx, constraints, jacobian_constr_T, jacobian_constr, constraint_start0, constraint_start1, hess_constr, hess_obj_fun, hess_lag_fun, grad_obj_fun, grad_lag_fun, w_inv, W_inv, Y, W_inv_Y, constraints_curr, grad_obj_curr, hess_lag_curr, jacobian_curr, grad_lag_curr, rho, gamma, sigma, y_inv, W_Y_inv, dir_x, dir_y, dir_w, b, stepsize_max, stepsize, merit_function, merit, grad_merit_fun, B, mu, mu_max, grad_merit_curr, merit_grad_dot_dir, dir_xw, infeasibility_tolerance, slack_tolerance, max_barrier_iterations, norm_grad_lag_start, norm_w_slack_start, iter_count; ASSERT( type( numVars, 'integer' ) ); ASSERT( nargs = 6 ); constraints := ConvertConstraints( constraintSet, variableBounding, numVars ); n := numVars; m := DimensionFast( constraints ); constraint_start0 := subs( {seq( x[i]=x_start[i], i=1..n )}, constraints ); constraint_start1 := map( x->max( 1, x ), constraint_start0 ); w_slack := Vector( map( evalf, constraint_start1 ), datatype=float[8] ); y_mult := Vector( m, .95, datatype=float[8] ); x_curr := Vector( map( evalf, x_start ), datatype=float[8] ); userinfo( 1, 'PrimalDualLogBarrier', `Minimize`, print(obj_function) ); userinfo( 1, 'PrimalDualLogBarrier', `subject to`, printConstraints(constraints, m) ); userinfo( 2, 'PrimalDualLogBarrier', `Starting point` ); userinfo( 2, 'PrimalDualLogBarrier', print(x_curr) ); userinfo( 3, 'PrimalDualLogBarrier', print(Y0), print(y_mult), print(W0), print(w_slack) ); grad_obj_fun := Vector( grad( obj_function, [seq( x[i], i=1..n )] ) ); hess_obj_fun := Matrix( hessian( obj_function, [seq( x[i], i=1..n )] ) ); jacobian_constr := Matrix( jacobian( constraints, [seq( x[i], i=1..n )] ) ); jacobian_constr_T := TransposeFast( jacobian_constr ); hess_constr := map( Matrix, map( hessian, constraints, [seq( x[i], i=1..n )] ) ); grad_lag_fun := VectorSubtractFast( grad_obj_fun, Vector( evalm( add( MultiplyVS( MatrixColumn( jacobian_constr_T, i ), k[i] ), i=1..m ) ) ) ); hess_lag_fun := MatrixSubtractFast( hess_obj_fun, Matrix( evalm( add( MultiplyMS( hess_constr[i], k[i] ), i=1..m ) ) ) ); merit := obj_function - mu*add( log( k[i] ), i=1..m ) + .5*B*add( (constraints[i]-k[i])^2, i=1..m ); merit_function := unapply( merit, seq( x[i], i=1..n ), seq( k[i], i=1..m ), mu, B ); grad_merit_fun := Vector( grad( merit, [seq( x[i], i=1..n ), seq( k[i], i=1..m )] ) ); userinfo( 3, 'PrimalDualLogBarrier', GRAD_F, print(TransposeFast( grad_obj_fun ) )); userinfo( 3, 'PrimalDualLogBarrier', HESSIAN_F, print(hess_obj_fun )); userinfo( 3, 'PrimalDualLogBarrier', JACOBIAN_OF_CONSTRAINTS, print(jacobian_constr )); userinfo( 3, 'PrimalDualLogBarrier', GRAD_LAGRANGIAN, print(TransposeFast( grad_lag_fun )) ); userinfo( 3, 'PrimalDualLogBarrier', HESSIAN_LAGRANGIAN, print(hess_lag_fun) ); userinfo( 3, 'PrimalDualLogBarrier', print(merit) ); userinfo( 3, 'PrimalDualLogBarrier', print(grad_merit_fun) ); B := .1; max_barrier_iterations := min( 200, 50 * n ); infeasibility_tolerance := .00001; slack_tolerance := .00001 * m; constraints_curr := evalf( EvaluateVector( constraints, x_curr, n ) ); grad_lag_curr := evalf( EvaluateDerivativeSum( grad_lag_fun, x_curr, y_mult, n, m ) ); rho := VectorSubtractFast( w_slack, constraints_curr ); norm_grad_lag_start := VectorNormFast( grad_lag_curr ); norm_w_slack_start := VectorNormFast( w_slack ); mu := InitialBarrierParameter( w_slack, constraints_curr, n, m ); mu_max := 100000 / m; iter_count := 1; to max_barrier_iterations while not ShouldStop( grad_lag_curr, norm_grad_lag_start, rho, norm_w_slack_start, w_slack, y_mult, infeasibility_tolerance, slack_tolerance ) do jacobian_curr := evalf( EvaluateDerivative( jacobian_constr, x_curr, n ) ); grad_obj_curr := evalf( EvaluateDerivative( grad_obj_fun, x_curr, n ) ); hess_lag_curr := evalf( EvaluateDerivativeSum( hess_lag_fun,x_curr, y_mult, n, m ) ); grad_merit_curr := evalf( EvaluateDerivativeSum(grad_merit_fun,x_curr,w_slack,n,m )); w_inv := map( x->1/x, w_slack ); y_inv := map( x->1/x, y_mult ); W_inv_Y := DiagMatrix( [seq( w_inv[i] * y_mult[i], i=1..m )], m ); W_Y_inv := DiagMatrix( [seq( w_slack[i] * y_inv[i], i=1..m )], m ); gamma := ScaledVectorAdd( w_inv, y_mult, mu, -1 ); b := VectorAddFast( MultiplyMV( W_inv_Y, rho ), gamma ); dir_xw := ComputeDirectionXW( grad_lag_curr, hess_lag_curr, jacobian_curr, W_inv_Y, b,grad_merit_curr, rho, constraints_curr, grad_obj_curr, n, convexityType ); dir_x := SubVectorFast( dir_xw, [1..n] ); dir_w := SubVectorFast( dir_xw, [n+1..n+m] ); dir_y := ComputeDirectionY( dir_w, W_inv_Y, b, rho ); merit_grad_dot_dir := evalf( DotProd( grad_merit_curr, dir_xw ) ); if merit_grad_dot_dir >= 0 then B := SetMeritPenaltyParameter( dir_x, grad_obj_curr, rho, mu, jacobian_curr, w_inv, m ); end if; stepsize := MeritStepSize( merit_function, n, m, x_curr, w_slack, y_mult, dir_x, dir_y, dir_w, dir_xw, merit_grad_dot_dir, mu, B ); if stepsize < .00000001 and iter_count >= 30 then printf("No solution found. Best approximation found (may not be feasible):"); return x_curr; elif stepsize < .000001 then B := 2 * B; end if; ScaledVectorAddInPlace( x_curr, dir_x, 1, stepsize ); ScaledVectorAddInPlace( y_mult, dir_y, 1, stepsize ) ; ScaledVectorAddInPlace( w_slack, dir_w, 1, stepsize ) ; userinfo( 3, 'PrimalDualLogBarrier', MU, mu, B_IS, B ); userinfo( 2, 'PrimalDualLogBarrier', print(TransposeFast( x_curr )) ); userinfo( 3, 'PrimalDualLogBarrier', print(dir_x), print(x_curr) ); userinfo( 3, 'PrimalDualLogBarrier', print(Y_DIR), print(dir_y), print(y_mult) ); userinfo( 3, 'PrimalDualLogBarrier', print(W_DIR), print(dir_w), print(w_slack) ); constraints_curr := evalf( EvaluateVector( constraints, x_curr, n ) ); grad_lag_curr := evalf( EvaluateDerivativeSum( grad_lag_fun, x_curr, y_mult, n, m )); rho := VectorSubtractFast( w_slack, constraints_curr ); userinfo( 4, 'PrimalDualLogBarrier', H( X ), constraints_curr, W, w_slack, INFEAS, rho ); mu := min( ComputeBarrierParameter( w_slack, y_mult, m, constraints_curr ), mu_max ); iter_count := iter_count + 1; end do; if convexityType = convex then userinfo( 2, 'PrimalDualLogBarrier', `Global optimum found at` ); else userinfo( 2, 'PrimalDualLogBarrier', `Local optimum found at` ); end if; { seq( x[i] = x_curr[i], i = 1..n) }; end proc: ShouldStop := proc(grad_lag_curr, norm_grad_lag_start, rho, norm_w_slack_start, w_slack, y_mult, infeas_tolerance, slack_tolerance ) local relative_primal_infeas, relative_dual_infeas, relative_comp_slack; relative_dual_infeas := VectorNormFast( grad_lag_curr ) / norm_grad_lag_start; relative_primal_infeas := VectorNormFast( rho ) / norm_w_slack_start; relative_comp_slack := DotProd( w_slack, y_mult ); userinfo( 3, 'PrimalDualLogBarrier', REL_DUAL_INFEAS, relative_dual_infeas, REL_PRIMAL_INFEAS, relative_primal_infeas, REL_COMP_SLACK, relative_comp_slack ); evalb( relative_dual_infeas <= infeas_tolerance and relative_primal_infeas <= infeas_tolerance and relative_comp_slack <= slack_tolerance ); end proc: InitialBarrierParameter := proc( w_slack, constraints_curr, n, m ) local big_M, mu_init; big_M := 200 / m; mu_init := big_M * max( VectorNormFast( w_slack ), VectorNormFast( constraints_curr ) ); userinfo( 2, 'Optimize', `Initial barrier parameter:`, mu_init ); mu_init; end proc: ComputeBarrierParameter := proc(w_slack, y_mult, m, constraints_curr) local dotprod_wy, mu_default, wy, min_wy, xi, i, scale; wy := seq( w_slack[i] * y_mult[i], i=1..m ); dotprod_wy := `+`( wy ); scale := .1; mu_default := scale * dotprod_wy / m; min_wy := min( wy ); xi := scale * min_wy / mu_default; if xi = 1 then mu_default; elif xi = 0 then 8 * mu_default; else min( 2, .95*( 1-xi )/xi )^3 * mu_default; end if; end proc: ComputeDirectionXW := proc(grad_lag_curr, hess_lag_curr, jacobian_curr, W_inv_Y, b, grad_merit_curr, rho, constraints_curr, grad_obj_curr, n, convexityType ) local W_inv_YA, A_trans_W_inv_YA, N, v, H, H_diag, lambda, dir_x, dir_w, i; W_inv_YA := MultiplyMM( W_inv_Y, jacobian_curr ); A_trans_W_inv_YA := MultiplyMM( TransposeFast( jacobian_curr ), W_inv_YA ); N := MatrixAddFast( hess_lag_curr, A_trans_W_inv_YA ); v := VectorSubtractFast( MultiplyMV( TransposeFast( jacobian_curr ), b ), grad_lag_curr ); if convexityType = nonconvex then H := hess_lag_curr; H_diag := DiagMatrix( [seq( max(abs( H[i, i] ),.01), i=1..n )], n ); lambda := 1.2; i := 1; to 20 while IsPosDefinite( N ) = false do ScaledMatrixAddInPlace( N, H_diag, 1, lambda ); userinfo(3,'Optimize', `Generalized Hessian adjusted by diagonal factor`, lambda ); if i < 10 then lambda := 2 * lambda; else lambda := 10 * lambda; end if; end do; end if; dir_x := LinearSolveFast( N, v ); dir_w := VectorSubtractFast( MultiplyMV( jacobian_curr, dir_x ), rho ); Vector( map( evalf, [dir_x, dir_w] ), datatype=float[8] ); end proc: ComputeDirectionY := proc(dir_w, W_inv_Y, b, rho ) local dir_y; dir_y := VectorSubtractFast( b, MultiplyMV( W_inv_Y, VectorAddFast( dir_w, rho ) ) ); Vector( map( evalf, dir_y ), datatype=float[8] ); end proc: ComputeDirectionW := proc(dir_y, W_Y_inv, gamma) #not currently called MultiplyMV( W_Y_inv, VectorSubtractFast( gamma, dir_y ) ); end proc: FeasibleStepSize := proc(dir_w, dir_y, w_slack, y_mult, m) local alpha_max, w_ratio, y_ratio, max_ratio, min_ratio; max_ratio := max( seq( -1*dir_w[i]/w_slack[i], i=1..m ), seq( -1*dir_y[i]/y_mult[i], i=1..m ) ); if max_ratio <= 1 then 0.95; else 0.95 / max_ratio; end if; end proc: MeritStepSize := proc(merit_fun, n, m, x_curr, w_slack, y_mult, dir_x, dir_y, dir_w, dir_xw, merit_grad_dot_dir, mu, B) local stepsize_max, stepsize_armijo, grad_merit_curr, xw_curr, direction_xw, grad_merit_xw, x_new, constraints_new; stepsize_max := FeasibleStepSize( dir_w, dir_y, w_slack, y_mult, m ); xw_curr := Vector( [x_curr, w_slack, mu, B] ); direction_xw := Vector( [dir_xw, 0, 0] ); stepsize_armijo := ConstrainedArmijoStepSize( merit_fun, n+m+2, xw_curr, direction_xw, merit_grad_dot_dir, stepsize_max ); userinfo( 3, 'NoName', "Max stepsize %g, Armijo stepsize %g\n", stepsize_max, stepsize_armijo ); stepsize_armijo; end proc: ConstrainedArmijoStepSize := proc(F, dim, point_curr, direction, grad_dot_direction, stepsize_max) local beta, sigma, m, point_new, improvement, threshold, beta_m; beta := .5; beta_m := stepsize_max; sigma := .001; m := 0; point_new := VectorAddFast( point_curr, direction ); improvement := F( seq( point_curr[i], i=1..dim ) ) - F( seq( point_new[i], i=1..dim ) ); threshold := -1 * sigma * grad_dot_direction; to 50 while improvement < threshold and abs( improvement ) > 0 do beta_m := beta * beta_m; point_new := ScaledVectorAdd( point_curr, direction, 1, beta_m ); improvement := F( seq( point_curr[i], i=1..dim ) ) - F( seq( point_new[i], i=1..dim ) ); threshold := evalf( beta * threshold ); m := m+1; end do; beta_m; end proc: ComplexArmijoStepSize := proc(F, dim, point_curr, direction, grad_dot_direction, stepsize_max) local beta, sigma, m, point_new, improvement, realImprovement, threshold, beta_m; beta := .5; beta_m := stepsize_max; sigma := .001; m := 0; point_new := VectorAddFast( point_curr, direction ); improvement := evalf( F( seq( point_curr[i], i=1..dim ) ) ) - evalf( F( seq( point_new[i], i=1..dim )) ); realImprovement := Re( improvement ); threshold := -1 * sigma * grad_dot_direction; to 50 while realImprovement < threshold and abs( realImprovement ) > 0 do beta_m := beta * beta_m; point_new := ScaledVectorAdd( point_curr, direction, 1, beta_m ); improvement := evalf( F( seq( point_curr[i], i=1..dim ) ) ) - evalf( F( seq( point_new[i], i=1..dim )) ); realImprovement := Re( improvement ); threshold := evalf( beta * threshold ); m := m+1; end do; beta_m; end proc: SetMeritPenaltyParameter := proc(dir_x, grad_obj_curr, rho, mu, jacobian_curr, w_inv, m) local v, k, jacob_T, squared_norm_rho; jacob_T := TransposeFast( jacobian_curr ); v := ScaledVectorAdd( grad_obj_curr, MultiplyMV( jacob_T, w_inv ), 1, -1*mu ); k := mu * DotProd( w_inv, rho ); squared_norm_rho := add( rho[i]^2, i=1..m ); 2 * ( DotProd( v, dir_x ) + k ) / squared_norm_rho; end proc: UnconstrainedNewton := proc( f, n, x_start, convexityType, datType ) local x_curr, direction, grad_fun, hess_fun, grad_start, H, H_diag, lambda, norm_grad_start, grad_curr, norm_grad_curr, hess_curr, hess_inv, neg_grad_curr, i, fx, f_curr, step_size, grad_dot_direction, max_step_size, max_iterations; if not (datType = 'float[8]' or datType = 'complex[8]') then error "datatype must be float[8] or complex[8] but received `%2`", 2, datType; end if; fx := unapply( f, seq( x[i], i=1..n ) ); grad_fun := Vector['row']( grad( f, [seq( x[i], i=1..n )] ) ); hess_fun := Matrix( hessian( f, [seq( x[i], i=1..n )] ) ); x_curr := Vector( map( evalf, x_start ), datatype=float[8] ); grad_start := evalf( map( x->Re(x), EvaluateDerivative( grad_fun, x_curr, n )) ); norm_grad_start := VectorNormFast( grad_start ); if norm_grad_start = 0 and convexityType = convex then return x_start; end if; userinfo( 2, 'UnconstrainedNewton', "Starting point" ); userinfo( 2, 'UnconstrainedNewton', print(x_curr) ); userinfo( 3, 'UnconstrainedNewton', Y0, y_mult, W0, w_slack ); userinfo( 3, 'UnconstrainedNewton', GRAD_F, grad_fun ); userinfo( 3, 'UnconstrainedNewton', HESSIAN_F, hess_fun ); grad_curr := grad_start; direction := MultiplyVS( grad_curr, -1 ); max_iterations := 100; if datType = 'float[8]' then to max_iterations while not ShouldTerminate( norm_grad_start, grad_curr, direction ) do neg_grad_curr := MultiplyVS( grad_curr, -1 ); hess_curr := evalf( EvaluateDerivative( hess_fun, x_curr, n ) ); H := hess_curr; if convexityType = nonconvex then H_diag := DiagMatrix( [seq( max(abs( H[i, i] ),.01), i=1..n )], n ); lambda := 1.2; i := 1; to 20 while IsPosDefinite( H ) = false do ScaledMatrixAddInPlace( H, H_diag, 1, lambda ); userinfo(3,'UnconstrainedNewton', `Hessian increased by diagonal factor`,lambda); if i < 10 then lambda := 2 * lambda; else lambda := 10 * lambda; end if; i := i + 1; end do; end if; direction := LinearSolveFast( H, neg_grad_curr ); grad_dot_direction := DotProd( grad_curr, direction ); max_step_size := .95; step_size := ConstrainedArmijoStepSize( fx, n, x_curr, direction, grad_dot_direction, max_step_size ); ScaledVectorAddInPlace( x_curr, direction, 1, step_size ); grad_curr := EvaluateDerivative( grad_fun, x_curr, n ); userinfo(2, 'UnconstrainedNewton', print(TransposeFast( x_curr )) ); userinfo(3, 'UnconstrainedNewton', X_DIR, print(direction), STEPSIZE, step_size, X, print(x_curr) ); end do: else # FOR COMPLEX MATRICES to max_iterations while not ShouldTerminate( norm_grad_start, grad_curr, direction ) do neg_grad_curr := evalf( MultiplyVS( map( x->Re(x), grad_curr), -1 ) ); hess_curr := evalf( map( x->Re(x), EvaluateDerivative( hess_fun, x_curr, n ) ) ); H := hess_curr; if convexityType = nonconvex then H_diag := DiagMatrix( [seq( max(abs( H[i, i] ),.01), i=1..n )], n ); lambda := 1.2; i := 1; to 20 while IsPosDefinite( H ) = false do ScaledMatrixAddInPlace( H, H_diag, 1, lambda ); userinfo(3,'UnconstrainedNewton', `Hessian increased by diagonal factor`, lambda); if i < 10 then lambda := 2 * lambda; else lambda := 10 * lambda; end if; i := i + 1; end do; end if; direction := LinearSolveFast( H, neg_grad_curr ); grad_dot_direction := DotProd( grad_curr, direction ); max_step_size := .95; step_size := ComplexArmijoStepSize( fx, n, x_curr, direction, grad_dot_direction, max_step_size ); ScaledVectorAddInPlace( x_curr, direction, 1, step_size ); grad_curr := evalf( EvaluateDerivative( grad_fun, x_curr, n ) ); userinfo(2, 'UnconstrainedNewton', print(TransposeFast( x_curr )) ); userinfo(3, 'UnconstrainedNewton', X_DIR, print(direction), STEPSIZE,step_size,X, print(x_curr) ); end do: end if; if convexityType = convex then userinfo( 2, 'UnconstrainedNewton', `Global optimum found at` ); else userinfo( 2, 'UnconstrainedNewton', `Local optimum found at` ); end if; { seq( x[i] = x_curr[i], i = 1..n) }; end proc: ShouldTerminate := proc( norm_grad_start, grad_curr, direction ) local norm_grad_curr, norm_direction, rel_grad_norm; norm_grad_curr := evalf( VectorNormFast( grad_curr ) ); norm_direction := evalf( VectorNormFast( direction ) ); rel_grad_norm := norm_grad_curr / norm_grad_start; userinfo( 3, 'UnconstrainedNewton', REL_GRAD_NORM, print(rel_grad_norm), REL_DIR_X, print(norm_direction) ); evalb( rel_grad_norm < .001 and norm_direction < .001 ); end proc: end module: # Saving the NonlinearProgramming Package into a Library # To save the NonlinearProgramming package into a library, remove the sharp character (#) from the libname assignment below and replace the example directory name "C:/mylib/nlp" with the name of the directory where you want the package saved. Then remove the # characters from the march and savelib commands. # #libname:="C:/mylib/nlp",libname; #march('create',libname[1],100); #savelib('NonlinearProgramming'); #with( NonlinearProgramming ); # # Now execute the worksheet by choosing from the menu bar # Edit->Execute->Worksheet