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_neqThis 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])
}
EODpp_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])
}
EODpp_addpm({At=>'Bot'},<<'EOD');
SEE ALSO
BUGS
Please report any bugs or suggestions at http://rt.cpan.org/
AUTHOR
Mark Grimes, <mgrimes@cpan.org>
COPYRIGHT AND LICENSE
This software is copyright (c) 2013 by Mark Grimes, <mgrimes@cpan.org>.
This is free software; you can redistribute it and/or modify it under the same terms as the Perl 5 programming language system itself.