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 = zeros(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 ); # Pars => 'dmat(m,m); dvec(m); # [o]sol(m); [o]lagr(q); [o]crval(); # amat(m,q); bvec(q); int meq(); # int [o]iact(q); int [o]nact(); # int [o]iter(s=2); [t] work(z); int [o]ierr(); my ( \$sol, \$lagr, \$crval, \$iact, \$nact, \$iter, \$ierr ) = qpgen2( \$Dmat->copy, \$dvec->copy, \$Amat->transpose->copy, \$avec->copy, \$meq, ); croak "qp: constraints are inconsistent, no solution!" if any(\$ierr == 1); croak "qp: matrix D in quadratic function is not positive definite!" if any(\$ierr == 2); croak "qp: some problem with mininization" if any(\$ierr); return { x => \$sol, lagr => \$lagr, crval => \$crval, iact => \$iact, nact => \$nact, iter => \$iter, ierr => \$ierr, }; } EOD pp_add_exported('', 'qp'); pp_addpm pp_line_numbers(__LINE__, q{ sub qp { my (\$Dmat, \$dvec, %args) = @_; my \$col = 0; my \$row = 1; my \$stack = 2; my \$n = pdl \$Dmat->dim(\$row); # D is an [n x n] matrix my \$s = pdl \$Dmat->dim(\$stack); # ... of \$s stacked problems # Default handling for A_eq and A_neq my \$A_eq = exists \$args{A_eq} ? \$args{A_eq} : zeroes(0, \$n); my \$A_neq = exists \$args{A_neq} ? \$args{A_neq} : zeroes(0, \$n); 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 # Default handling for a_eq and a_neq # These have to be [?x1] matrixes to ensure threading works if glue, # otherwise they are relicated not just attached. my \$a_eq = exists \$args{a_eq} ? \$args{a_eq} : zeroes(1, \$m); my \$a_neq = exists \$args{a_neq} ? \$args{a_neq} : zeroes(1, \$p); check_dimmensions( 'Dmat', \$Dmat, \$n, \$n, \$s ); check_dimmensions( 'dvec', \$dvec, \$n, \$s ); check_dimmensions( 'A_eq', \$A_eq, \$m, \$n, \$s ); check_dimmensions( 'a_eq', \$a_eq, \$m, \$s ); check_dimmensions( 'A_neq', \$A_neq, \$p, \$n, \$s ); if( \$p == 0 ){ # If a_neq has zero elems (per stack), then it will be [0x1] check_dimmensions( 'a_neq', \$a_neq, 1, \$p, \$s ); } else { # otherwise it will be a [p] vector check_dimmensions( 'a_neq', \$a_neq, \$p, \$s ); } # Combine _eq and _neq, use meq to specifiy number of equality constratins my \$A = \$A_eq->glue( 0, \$A_neq ); my \$a = \$a_eq->glue( 0, \$a_neq ); my \$meq = \$A_eq->dim(\$col); # Pars => 'dmat(m,m); dvec(m); # [o]sol(m); [o]lagr(q); [o]crval(); # amat(m,q); bvec(q); int meq(); # int [o]iact(q); int [o]nact(); # int [o]iter(s=2); [t] work(z); int [io]ierr(); my ( \$sol, \$lagr, \$crval, \$iact, \$nact, \$iter, \$ierr ) = qpgen2( \$Dmat->copy, \$dvec->copy, \$A->transpose->copy, \$a->copy, \$meq, ); croak "qp: constraints are inconsistent, no solution!" if any(\$ierr == 1); croak "qp: matrix D in quadratic function is not positive definite!" if any(\$ierr == 2); croak "qp: some problem with mininization" if any(\$ierr); return { x => \$sol, lagr => \$lagr, crval => \$crval, iact => \$iact, nact => \$nact, iter => \$iter, ierr => \$ierr, }; } sub check_dimmensions { my (\$name, \$pdl, @dims) = @_; pop @dims if \$dims[-1] == 1; # remove last dim if a dimension of 1 my \$got = pdl \$pdl->dims; my \$expected = pdl @dims; croak( sprintf( "dimmension check failed for %s (is %s)\nexpected [%s]\n%s", \$name, \$pdl->info, join(',',@dims), qp_usage() ) ) unless all \$got == \$expected; } sub qp_usage { return q{ usage: qp( Dmat [n x n], dvec [n], A_eq => [n x m], a_eq => [m], a_neq => [n x p], a_neq => [p] ) where A_eq and a_eq are the equality constraints and A_neq and a_neq are the inequality constraints (>=) All inputs can add one addition dimmension to "stack" multiple Optimization problems. The returned solution (x) will have stack dimmensions. }; } }); 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.