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.