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.