GALAHAD DPS package#

purpose#

The dps package constructs a symmetric, positive definite matrix $M$ from a given $H$ so that $H$ is is diagonal in the norm $\|v\|_M = \sqrt{v^T M v}$ induced by $M$, and consequently minimizers of trust-region and regularized quadratic subproblems may be computed efficiently. The aim is either to minimize the quadratic objective function

\[q(x) = f + g^T x + \frac{1}{2} x^T H x,\]

where the vector $x$ is required to satisfy the ellipsoidal trust-region constraint $\|x\|_{M} \leq \Delta$, or to minimize the regularized quadratic objective

\[r(x) = q(x) + \frac{\sigma}{p} \|x\|_M^p,\]

where the radius $\Delta > 0$, the weight $\sigma > 0$, and the power $p \geq 2$. A factorization of the matrix $H$ will be required, so this package is most suited for the case where such a factorization, either dense or sparse, may be found efficiently.

See Section 4 of $GALAHAD/doc/dps.pdf for additional details.

method#

The required solution $x_*$ necessarily satisfies the optimality condition $H x_* + \lambda_* M x_* + g = 0$, where $\lambda_* \geq 0$ is a Lagrange multiplier that corresponds to the constraint $\|x\|_M \leq \Delta$ in the trust-region case, and is given by $\lambda_* = \sigma \|x_*\|^{p-2}$ for the regularization problem involve $r(x)$. In addition $H + \lambda_* M$ will be positive semi-definite; in most instances it will actually be positive definite, but in special “hard” cases singularity is a possibility.

The matrix $H$ is decomposed as

\[H = P L D L^T P^T\]

by calling the matrix factorization package SLS. Here $P$ is a permutation matrix, $L$ is unit lower triangular and $D$ is block diagonal, with blocks of dimension at most two. The spectral decomposition of each diagonal block of $D$ is computed, and each eigenvalue $\theta$ is replaced by $\max ( | \theta | , \theta_{\min} ) $, where $\theta_{\min}$ is a positive user-supplied value. The resulting block diagonal matrix is $B$, from which we define the modified-absolute-value

\[M = P L B L^T P^T;\]

an alternative due to Goldfarb uses instead the simpler

\[M = P L L^T P^T.\]

Given the factors of $H$ (and $M$), the required solution is found by making the change of variables $y = B^{1/2} L^T P^T x$ (or $y = L^T P^T x$ in the Goldfarb case) which results in ``diagonal’’ trust-region and regularization subproblems, whose solution may be easily obtained suing a Newton or higher-order iteration of a resulting “secular” equation. If subsequent problems, for which $H$ and $g$ are unchanged, are to be attempted, the existing factorization and solution may easily be exploited.

The dominant cost is that for the factorization of the symmetric, but potentially indefinite, matrix $H$ using the package SLS.

references#

The method is described in detail for the trust-region case in

N. I. M. Gould and J. Nocedal, ``The modified absolute-value factorization for trust-region minimization’’. In ``High Performance Algorithms and Software in Nonlinear Optimization’’ (R. De Leone, A. Murli, P. M. Pardalos and G. Toraldo, eds.), Kluwer Academic Publishers, (1998) 225–241,

while the adaptation for the regularization case is obvious. The method used to solve the diagonal trust-region and regularization subproblems are as given by

H. S. Dollar, N. I. M. Gould and D. P. Robinson, ``On solving trust-region and other regularised subproblems in optimization’’. Mathematical Programming Computation 2(1) (2010) 21–57

with simplifications due to the diagonal Hessian.

introduction to function calls#

To solve a given problem, functions from the dps package must be called in the following order:

dps_initialize - provide default control parameters and set up initial data structures
dps_read_specfile (optional) - override control values by reading replacement values from a file
dps_import - import control and matrix data structures
dps_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
one of
- dps_solve_tr_problem - solve the trust-region problem (1)
- dps_solve_rq_problem - solve the regularized-quadratic problem (2)
optionally one of
- dps_resolve_tr_problem - resolve the trust-region problem (1) when the non-matrix data has changed
- dps_resolve_rq_problem - resolve the regularized-quadratic problem (2) when the non-matrix data has changed
dps_information (optional) - recover information about the solution and solution process
dps_terminate - deallocate data structures

See the examples section for illustrations of use.

parametric real type T#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision), Float64 (double precision) or, if supported, Float128 (quadruple precision).

callable functions#

    function dps_initialize(T, data, control, status)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a structure containing control information (see dps_control_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

0
The initialization was successful.

    function dps_read_specfile(T, control, specfile)

Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/dps/DPS.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/dps.pdf for a list of how these keywords relate to the components of the control structure.

Parameters:

control	is a structure containing control information (see dps_control_type)
specfile	is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file

    function dps_import(T, control, data, status, n,
                        H_type, ne, H_row, H_col, H_ptr)

Import problem data into internal storage prior to solution.

Parameters:

control	is a structure whose members provide control parameters for the remaining procedures (see dps_control_type)
data	holds private internal data
status	is a scalar variable of type Int32 that gives the exit status from the package. Possible values are: 1 The import was successful, and the package is ready for the solve phase -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated.
n	is a scalar variable of type Int32 that holds the number of variables
H_type	is a one-dimensional array of type Vararg{Cchar} that specifies the symmetric storage scheme used for the Hessian. It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense’; lower or upper case variants are allowed
ne	is a scalar variable of type Int32 that holds the number of entries in the lower triangular part of H in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes.
H_row	is a one-dimensional array of size ne and type Int32 that holds the row indices of the lower triangular part of H in the sparse co-ordinate storage scheme. It need not be set for any of the other three schemes, and in this case can be C_NULL
H_col	is a one-dimensional array of size ne and type Int32 that holds the column indices of the lower triangular part of H in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense or diagonal storage schemes are used, and in this case can be C_NULL
H_ptr	is a one-dimensional array of size n+1 and type Int32 that holds the starting position of each row of the lower triangular part of H, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be C_NULL

    function dps_reset_control(T, control, data, status)

Reset control parameters after import if required.

Parameters:

control

is a structure whose members provide control parameters for the remaining procedures (see dps_control_type)

data

holds private internal data

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are:

1
The import was successful, and the package is ready for the solve phase

    function dps_solve_tr_problem(T, data, status, n, ne, H_val, c, f,
                                  radius, x)

Find the global minimizer of the trust-region problem (1).

Parameters:

data	holds private internal data
status	is a scalar variable of type Int32 that gives the exit status from the package. Possible values are: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -40 An error has occured when building the preconditioner.
n	is a scalar variable of type Int32 that holds the number of variables
ne	is a scalar variable of type Int32 that holds the number of entries in the lower triangular part of the Hessian matrix $H$.
H_val	is a one-dimensional array of size ne and type T that holds the values of the entries of the lower triangular part of the Hessian matrix $H$ in any of the available storage schemes.
c	is a one-dimensional array of size n and type T that holds the linear term $c$ in the objective function. The j-th component of `c`, j = 1, … , m, contains $c_j$.
f	is a scalar variable pointer of type T that holds the value of the holds the constant term $f$ in the objective function.
radius	is a scalar variable pointer of type T that holds the value of the trust-region radius, $\Delta > 0$.
x	is a one-dimensional array of size n and type T that holds the values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.

    function dps_solve_rq_problem(T, data, status, n, ne, H_val, c, f,
                                  power, weight, x)

Find the global minimizer of the regularized-quadartic problem (2).

Parameters:

data	holds private internal data
status	is a scalar variable of type Int32 that gives the exit status from the package. Possible values are: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -40 An error has occured when building the preconditioner.
n	is a scalar variable of type Int32 that holds the number of variables
ne	is a scalar variable of type Int32 that holds the number of entries in the lower triangular part of the Hessian matrix $H$.
H_val	is a one-dimensional array of size ne and type T that holds the values of the entries of the lower triangular part of the Hessian matrix $H$ in any of the available storage schemes.
c	is a one-dimensional array of size n and type T that holds the linear term $c$ in the objective function. The j-th component of `c`, j = 1, … , m, contains $c_j$.
f	is a scalar variable pointer of type T that holds the value of the holds the constant term $f$ in the objective function.
weight	is a scalar variable pointer of type T that holds the value of the regularization weight, $\sigma > 0$.
power	is a scalar variable pointer of type T that holds the value of the regularization power, $p \geq 2$.
x	is a one-dimensional array of size n and type T that holds the values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.

    function dps_resolve_tr_problem(T, data, status, n, c, f, radius, x)

Find the global minimizer of the trust-region problem (1) if some non-matrix components have changed since a call to dps_solve_tr_problem.

Parameters:

data	holds private internal data
status	is a scalar variable of type Int32 that gives the exit status from the package. Possible values are: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -16 The problem is so ill-conditioned that further progress is impossible.
n	is a scalar variable of type Int32 that holds the number of variables
c	is a one-dimensional array of size n and type T that holds the linear term $c$ in the objective function. The j-th component of `c`, j = 1, … , m, contains $c_j$.
f	is a scalar variable pointer of type T that holds the value of the constant term $f$ in the objective function.
radius	is a scalar variable pointer of type T that holds the value of the trust-region radius, $\Delta > 0$.
x	is a one-dimensional array of size n and type T that holds the values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.

    function dps_resolve_rq_problem(T, data, status, n, c, f, power, weight, x)

Find the global minimizer of the regularized-quadartic problem (2) if some non-matrix components have changed since a call to dps_solve_rq_problem.

Parameters:

data	holds private internal data
status	is a scalar variable of type Int32 that gives the exit status from the package. Possible values are: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -16 The problem is so ill-conditioned that further progress is impossible.
n	is a scalar variable of type Int32 that holds the number of variables
c	is a one-dimensional array of size n and type T that holds the linear term $c$ in the objective function. The j-th component of `c`, j = 1, … , m, contains $c_j$.
f	is a scalar variable pointer of type T that holds the value of the holds the constant term $f$ in the objective function.
weight	is a scalar variable pointer of type T that holds the value of the regularization weight, $\sigma > 0$.
power	is a scalar variable pointer of type T that holds the value of the regularization power, $p \geq 2$.
x	is a one-dimensional array of size n and type T that holds the values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.

    function dps_information(T, data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see dps_inform_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

0
The values were recorded successfully

    function dps_terminate(T, data, control, inform)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a structure containing control information (see dps_control_type)
inform	is a structure containing output information (see dps_inform_type)

available structures#

dps_control_type structure#

    struct dps_control_type{T}
      f_indexing::Bool
      error::Int32
      out::Int32
      problem::Int32
      print_level::Int32
      new_h::Int32
      taylor_max_degree::Int32
      eigen_min::T
      lower::T
      upper::T
      stop_normal::T
      stop_absolute_normal::T
      goldfarb::Bool
      space_critical::Bool
      deallocate_error_fatal::Bool
      problem_file::NTuple{31,Cchar}
      symmetric_linear_solver::NTuple{31,Cchar}
      prefix::NTuple{31,Cchar}
      sls_control::sls_control_type{T}

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

Int32 error

unit for error messages

Int32 out

unit for monitor output

Int32 problem

unit to write problem data into file problem_file

Int32 print_level

controls level of diagnostic output

Int32 new_h

how much of $H$ has changed since the previous call. Possible values are

0 unchanged
1 values but not indices have changed
2 values and indices have changed

Int32 taylor_max_degree

maximum degree of Taylor approximant allowed

T eigen_min

smallest allowable value of an eigenvalue of the block diagonal factor of $H$

T lower

lower and upper bounds on the multiplier, if known

T upper

see lower

T stop_normal

stop trust-region solution when $| ||x||_M - \delta | \leq$ max( .stop_normal \* delta, .stop_absolute_normal )

T stop_absolute_normal

see stop_normal

Bool goldfarb

use the Goldfarb variant of the trust-region/regularization norm rather than the modified absolute-value version

Bool space_critical

if space is critical, ensure allocated arrays are no bigger than needed

Bool deallocate_error_fatal

exit if any deallocation fails

char problem_file[31]

name of file into which to write problem data

char symmetric_linear_solver[31]

symmetric (indefinite) linear equation solver

NTuple{31,Cchar} prefix

all output lines will be prefixed by prefix(2:LEN(TRIM(.prefix))-1) where prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’

struct sls_control_type sls_control

control parameters for the Cholesky factorization and solution

dps_time_type structure#

    struct dps_time_type{T}
      total::T
      analyse::T
      factorize::T
      solve::T
      clock_total::T
      clock_analyse::T
      clock_factorize::T
      clock_solve::T

detailed documentation#

time derived type as a Julia structure

components#

T total

total CPU time spent in the package

T analyse

CPU time spent reordering H prior to factorization.

T factorize

CPU time spent factorizing H.

T solve

CPU time spent solving the diagonal model system.

T clock_total

total clock time spent in the package

T clock_analyse

clock time spent reordering H prior to factorization

T clock_factorize

clock time spent factorizing H

T clock_solve

clock time spent solving the diagonal model system

dps_inform_type structure#

    struct dps_inform_type{T}
      status::Int32
      alloc_status::Int32
      mod_1by1::Int32
      mod_2by2::Int32
      obj::T
      obj_regularized::T
      x_norm::T
      multiplier::T
      pole::T
      hard_case::Bool
      bad_alloc::NTuple{81,Cchar}
      time::dps_time_type{T}
      sls_inform::sls_inform_type{T}

detailed documentation#

inform derived type as a Julia structure

components#

Int32 status

return status. See DPS_solve for details

Int32 alloc_status

STAT value after allocate failure.

Int32 mod_1by1

the number of 1 by 1 blocks from the factorization of H that were modified when constructing $M$

Int32 mod_2by2

the number of 2 by 2 blocks from the factorization of H that were modified when constructing $M$

T obj

the value of the quadratic function

T obj_regularized

the value of the regularized quadratic function

T x_norm

the M-norm of the solution

T multiplier

the Lagrange multiplier associated with the constraint/regularization

T pole

a lower bound max(0,-lambda_1), where lambda_1 is the left-most eigenvalue of $(H,M)$

Bool hard_case

has the hard case occurred?

NTuple{81,Cchar} bad_alloc

name of array that provoked an allocate failure

struct dps_time_type time

time information

struct sls_inform_type sls_inform

information from SLS

Index

example calls#

This is an example of how to use the package to solve a trust-region subproblem; the code is available in $GALAHAD/src/dps/Julia/test_dps.jl . A variety of supported Hessian and constraint matrix storage formats are shown.

# test_dps.jl
# Simple code to test the Julia interface to DPS

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_dps(::Type{T}) where T
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{dps_control_type{T}}()
  inform = Ref{dps_inform_type{T}}()

  # Set problem data
  n = 3 # dimension of H
  m = 1 # dimension of A
  H_ne = 4 # number of elements of H
  H_dense_ne = 6 # number of elements of H
  H_row = Cint[1, 2, 3, 3]  # row indices, NB lower triangle
  H_col = Cint[1, 2, 3, 1]
  H_ptr = Cint[1, 2, 3, 5]
  H_val = T[1.0, 2.0, 3.0, 4.0]
  H_dense = T[1.0, 0.0, 2.0, 4.0, 0.0, 3.0]
  f = T(0.96)
  radius = one(T)
  half_radius = T(0.5)
  c = T[0.0, 2.0, 0.0]

  st = ' '
  status = Ref{Cint}()
  x = zeros(T, n)

  @printf(" Fortran sparse matrix indexing\n\n")
  @printf(" basic tests of storage formats\n\n")

  for storage_type in 1:3
    # Initialize DPS
    dps_initialize(T, data, control, status)

    # Set user-defined control options
    @reset control[].f_indexing = true # fortran sparse matrix indexing
    @reset control[].symmetric_linear_solver = galahad_linear_solver("sytr")

    # sparse co-ordinate storage
    if storage_type == 1
      st = 'C'
      # import the control parameters and structural data
      dps_import(T, control, data, status, n,
                 "coordinate", H_ne, H_row, H_col, C_NULL)

      # solve the problem
      dps_solve_tr_problem(T, data, status, n, H_ne, H_val,
                           c, f, radius, x)
    end

    # sparse by rows
    if storage_type == 2
      st = 'R'
      # import the control parameters and structural data
      dps_import(T, control, data, status, n,
                 "sparse_by_rows", H_ne, C_NULL, H_col, H_ptr)

      dps_solve_tr_problem(T, data, status, n, H_ne, H_val,
                           c, f, radius, x)
    end

    # dense
    if storage_type == 3
      st = 'D'
      # import the control parameters and structural data
      dps_import(T, control, data, status, n,
                 "dense", H_ne, C_NULL, C_NULL, C_NULL)

      dps_solve_tr_problem(T, data, status, n, H_dense_ne, H_dense,
                           c, f, radius, x)
    end

    dps_information(T, data, inform, status)
    @printf("format %c: DPS_solve_problem exit status   = %1i, f = %.2f\n",
            st, inform[].status, inform[].obj)

    # sparse co-ordinate storage
    if storage_type == 1
      st = 'C'
      # solve the problem
      dps_resolve_tr_problem(T, data, status, n,
                             c, f, half_radius, x)
    end

    # sparse by rows
    if storage_type == 2
      st = 'R'
      dps_resolve_tr_problem(T, data, status, n,
                             c, f, half_radius, x)
    end

    # dense
    if storage_type == 3
      st = 'D'
      dps_resolve_tr_problem(T, data, status, n,
                             c, f, half_radius, x)
    end

    dps_information(T, data, inform, status)
    @printf("format %c: DPS_resolve_problem exit status = %1i, f = %.2f\n",
            st, inform[].status, inform[].obj)

    # @printf("x: ")
    # for i = 1:n+m
    #   @printf("%f ", x[i])
    # end

    # Delete internal workspace
    dps_terminate(T, data, control, inform)
  end

  return 0
end

@testset "DPS" begin
  @test test_dps(Float32) == 0
  @test test_dps(Float64) == 0
  @test test_dps(Float128) == 0
end