NAME PDL::Opt::QP - Quadratic programming solver for PDL SYNOPSIS use PDL::Opt::QP; ... DESCRIPTION ... 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'); pp_addpm({At=>'Bot'},<<'EOD'); sub qp { 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); } die("Dmat is not square!") if $n != $Dmat->dim(0); # Check D is [n x n] die("Dmat and dvec are incompatible!") if $n != $dvec->nelem; # Check d is [n] die("Amat and dvec are incompatible!") if $n != $Amat->dim(1); # Check A is [n x _] die("Amat and bvec are incompatible!") if $q != $bvec->nelem; # Check A is [_ x q] die("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_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.