Mark Grimes > PDL-Opt-QP-0.27 > PDL::Opt::QP

Download:
PDL-Opt-QP-0.27.tar.gz

Dependencies

Annotate this POD

View/Report Bugs
Module Version: 0.27   Source  

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 http://rt.cpan.org/

AUTHOR ^

Mark Grimes, <mgrimes@cpan.org>

COPYRIGHT AND LICENSE ^

This software is copyright (c) 2014 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.

syntax highlighting: