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: A' x = a B x >= b This routine solves the quadratic programming optimization problem minimize f(x) = 0.5 x' D x - d' x x optionally constrained by: A' x = a B x >= b .... more docs to come .... }); pp_add_exported('', 'qp_orig'); pp_addpm({At=>'Bot'},<<'EOD'); sub qp_orig { my (\$Dmat, \$dvec, \$Amat, \$bvec, \$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( ! defined \$bvec ){ # if b is undef, create it \$bvec = zeros(\$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 bvec are incompatible!") if \$q != \$bvec->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 my \$res = qpgen2( \$Dmat->copy, \$dvec->copy, \$n, \$n, \$sol, \$lagr, \$crval, \$Amat->transpose->copy, \$bvec->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, \$Amat, \$avec, \$Bmat, \$bvec) = @_; my \$n = pdl \$Dmat->dim(1); # D is an [n x n] matrix my \$q = pdl \$Amat->dim(0) + \$Bmat->dim(0); # A is an [n x q] matrix if( ! defined \$bvec ){ # if b is undef, create it \$bvec = zeros(\$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("Bmat and dvec are incompatible!") if \$n != \$Bmat->dim(1); # Check B is [n x _] croak("Amat and avec are incompatible!") if \$Amat->dim(0) != \$avec->nelem; # Check A is [_ x q] croak("Bmat and bvec are incompatible!") if \$Bmat->dim(0) != \$bvec->nelem; # Check B is [_ x q] 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 = \$Amat->glue( 0, \$Bmat )->copy; my \$a = \$avec->glue( 0, \$bvec )->copy; # ->flat? my \$meq = \$Amat->dim(0); print "\\$A' = ", \$A->transpose; print "\\$a = ", \$a; print "n = \$n"; print "q = \$q"; print "meq = \$meq"; print "\n"; my \$res = qpgen2( \$Dmat->copy, \$dvec->copy, \$n, \$n, \$sol, \$lagr, \$crval, \$A->transpose, \$a, \$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.