NAME PDL::Opt::QP - Quadratic programming solver for PDL SYNOPSIS use PDL; use PDL::NiceSlice; use PDL::Opt::QP; my \$mu = pdl(q[ 0.0427 0.0015 0.0285 ])->transpose; # [ n x 1 ] my \$mu_0 = 0.0427; my \$dmat = pdl q[ 0.0100 0.0018 0.0011 ; 0.0018 0.0109 0.0026 ; 0.0011 0.0026 0.0199 ]; my \$dvec = zeros(3); my \$amat = \$mu->glue( 0, ones( 1, 3 ) )->copy; my \$bvec = pdl(\$mu_0)->glue( 1, ones(1) )->flat; my \$meq = pdl(2); my \$sol = qp( \$dmat, \$dvec, \$amat, \$bvec, \$meq ); say "Solution: ", \$sol->{x}; DESCRIPTION This routine uses Goldfarb/Idnani algorithm to solve the following minimization problem: minimize f(x) = 0.5 * x' D x - d' x x optionally constrained by: Aeq' x = a_eq Aneq x >= b_neq This routine solves the quadratic programming optimization problem minimize f(x) = 0.5 x' D x - d' x x optionally constrained by: Aeq' x = a_eq Aneq x >= b_neq .... more docs to come .... }); pp_add_exported('', 'qp_orig'); pp_addpm({At=>'Bot'},<<'EOD'); sub qp_orig { my (\$Dmat, \$dvec, \$Amat, \$avec, \$meq) = @_; my \$n = pdl \$Dmat->dim(1); # D is an [n x n] matrix my \$q = pdl \$Amat->dim(0); # A is an [n x q] matrix if( \$avec->isnull ){ \$avec = pdl()->reshape(1,\$q); } croak("Dmat is not square!") if \$n != \$Dmat->dim(0); # Check D is [n x n] croak("Dmat and dvec are incompatible!") if \$n != \$dvec->nelem; # Check d is [n] croak("Amat and dvec are incompatible!") if \$n != \$Amat->dim(1); # Check A is [n x _] croak("Amat and avec are incompatible!") if \$q != \$avec->nelem; # Check A is [_ x q] croak("Value of meq is invalid!") if (\$meq > \$q) || (\$meq < 0 ); my \$iact = zeros(\$q); # Store which constraints are active my \$nact = pdl(0); # Store number of active constraints my \$r = \$n < \$q ? \$n : \$q; # Used to size work space my \$sol = zeros(\$n->sclr); # Store the solution [n] my \$lagr = zeros(\$q->sclr); # Store the Lagranges for the constraints my \$crval= pdl(0); # Value at min my \$work = zeros((2*\$n+\$r*(\$r+5)/2+2*\$q+1)->sclr); # Work space my \$iter = zeros(2); # Store info about interations my \$ierr = pdl(0); # Input: 1=Factorized; Output: error flag # print "\\$A = ", \$Amat; # print "\\$a = \$avec\n"; # print "n = \$n\n"; # print "q = \$q\n"; # print "meq = \$meq\n"; my \$res = qpgen2( \$Dmat->copy, \$dvec->copy, \$n, \$n, \$sol, \$lagr, \$crval, \$Amat->transpose->copy, \$avec->copy, \$n, \$q, \$meq, \$iact, \$nact, \$iter, \$work, \$ierr ); croak "qp: constraints are inconsistent, no solution!" if \$ierr->sclr == 1; croak "qp: matrix D in quadratic function is not positive definite!" if \$ierr->sclr == 2; croak "qp: some problem with mininization" if \$ierr->sclr; return { x => \$sol, lagr => \$lagr, crval => \$crval, iact => \$iact, nact => \$nact, iter => \$iter, ierr => \$ierr, }; # TODO: process/return the results # # From R implementation: # # list(solution=res1\$sol, # value=res1\$crval, # unconstrained.solution=res1\$dvec, # iterations=res1\$iter, # Lagrangian = res1\$lagr, # iact=res1\$iact[1:res1\$nact]) } EOD pp_add_exported('', 'qp'); pp_addpm({At=>'Bot'},<<'EOD'); sub qp { my (\$Dmat, \$dvec, \$A_eq, \$a_eq, \$A_neq, \$a_neq) = @_; my \$col = 0; my \$row = 1; my \$n = pdl \$Dmat->dim(\$row); # D is an [n x n] matrix my \$m = pdl \$A_eq->dim(\$col); # A is an [n x m] matrix my \$p = pdl \$A_neq->dim(\$col); # A is an [n x p] matrix my \$q = \$m + \$p; # A is an [n x q] matrix if( \$A_eq->isnull ){ \$A_eq = pdl()->reshape(0,\$n); } if( \$a_eq->isnull ){ \$a_eq = pdl()->reshape(1,\$m); } if( \$A_neq->isnull ){ \$A_neq = pdl()->reshape(0,\$n); } if( \$a_neq->isnull ){ \$a_neq = pdl()->reshape(1,\$p); } croak("dimmension check failed: Dmat [n x n*] is not square") if \$Dmat->dim(\$col) != \$n; croak("dimmension check failed: Dmat [n x n] and dvec [n* x 1]") if \$dvec->nelem != \$n; croak("dimmension check failed: A_eq [n* x m] and a_eq [n x 1]") if \$A_eq->dim(\$row) != \$n; croak("dimmension check failed: A_eq [n x m] and a_eq [m* x 1]") if \$a_eq->nelem != \$m; croak("dimmension check failed: A_neq [n* x p] and a_neq [p x 1]") if \$A_neq->dim(\$row) != \$n; croak("dimmension check failed: A_neq [n x p] and a_neq [p* x 1]") if \$a_neq->nelem != \$p; my \$iact = zeros(\$q); # Store which constraints are active my \$nact = pdl(0); # Store number of active constraints my \$r = \$n < \$q ? \$n : \$q; # Used to size work space my \$sol = zeros(\$n->sclr); # Store the solution [n] my \$lagr = zeros(\$q->sclr); # Store the Lagranges for the constraints my \$crval= pdl(0); # Value at min my \$work = zeros((2*\$n+\$r*(\$r+5)/2+2*\$q+1)->sclr); # Work space my \$iter = zeros(2); # Store info about interations my \$ierr = pdl(0); # Input: 1=Factorized; Output: error flag my \$A = \$A_eq->glue( 0, \$A_neq )->copy; my \$a = \$a_eq->glue( 0, \$a_neq )->copy; # ->flat? my \$meq = \$A_eq->dim(0); # print "\\$A = ", \$A; # print "\\$a = \$a\n"; # print "n = \$n\n"; # print "q = \$q\n"; # print "meq = \$meq\n"; my \$res = qpgen2( \$Dmat->copy, \$dvec->copy, \$n, \$n, \$sol, \$lagr, \$crval, \$A->transpose->copy, \$a->copy, \$n, \$q, \$meq, \$iact, \$nact, \$iter, \$work, \$ierr ); croak "qp: constraints are inconsistent, no solution!" if \$ierr->sclr == 1; croak "qp: matrix D in quadratic function is not positive definite!" if \$ierr->sclr == 2; croak "qp: some problem with mininization" if \$ierr->sclr; return { x => \$sol, lagr => \$lagr, crval => \$crval, iact => \$iact, nact => \$nact, iter => \$iter, ierr => \$ierr, }; # TODO: process/return the results # # From R implementation: # # list(solution=res1\$sol, # value=res1\$crval, # unconstrained.solution=res1\$dvec, # iterations=res1\$iter, # Lagrangian = res1\$lagr, # iact=res1\$iact[1:res1\$nact]) } EOD pp_addpm({At=>'Bot'},<<'EOD'); SEE ALSO PDL, PDL::Opt::NonLinear BUGS Please report any bugs or suggestions at AUTHOR Mark Grimes, COPYRIGHT AND LICENSE This software is copyright (c) 2013 by Mark Grimes, . This is free software; you can redistribute it and/or modify it under the same terms as the Perl 5 programming language system itself.