GALAHAD GLTR package#

purpose#

The gltr package uses a Krylov-subspace iteration to find an approximation of the global minimizer of quadratic objective function within an ellipsoidal region; this is commonly known as the trust-region subproblem. The aim is 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$, where the $M$-norm of $x$ is defined to be $\|x\|_{M} = \sqrt{x^T M x}$, and where the radius $\Delta > 0$. The method may be suitable for large problems as no factorization of $H$ is required. Reverse communication is used to obtain matrix-vector products of the form $H z$ and $M^{-1} z.$

See Section 4 of $GALAHAD/doc/gltr.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 corresponding to the constraint $\|x\|_M \leq \Delta$. In addition, the matrix $H + \lambda M$ will be positive definite.

The method is iterative. Starting with the vector $M^{-1} g$, a matrix of Lanczos vectors is built one column at a time so that the $k$-th column is generated during iteration $k$. These columns span a so-called Krylov space. The resulting $n$ by $k$ matrix $Q_k $ has the property that $Q_{k}^T H Q_k^{} = T_{k}^{}$, where $T_k$ is tridiagonal. An approximation to the required solution may then be expressed formally as

\[x_{k+1} = Q_k y_k\]

where $y_k $ solves the ``tridiagonal’’ subproblem of minimizing

\[\frac{1}{2} y^T T_k y + \|g\|_{M^{-1} } e_{1}^T y \;\;\; \mbox{subject to the constraint}\;\; \|y\|_2 \leq \Delta,\;\;\mbox{(1)}\]

and where $e_1$ is the first unit vector.

If the solution to (1) lies interior to the constraint, the required solution $x_{k+1}$ may simply be found as the $k$-th (preconditioned) conjugate-gradient iterate. This solution can be obtained without the need to access the whole matrix $Q_k$. These conjugate-gradient iterates increase in $M$-norm, and thus once one of them exceeds $\Delta$ in $M$-norm, the solution must occur on the constraint boundary. Thereafter, the solution to (1) is less easy to obtain, but an efficient inner iteration to solve (1) is nonetheless achievable because $T_k $ is tridiagonal. It is possible to observe the optimality measure $\|H x + \lambda M x + g\|_{M^{-1}}$ without computing $x_{k+1}$, and thus without needing $Q_k $. Once this measure is sufficiently small, a second pass is required to obtain the estimate $x_{k+1} $ from $y_k $. As this second pass is an additional expense, a record is kept of the optimal objective function values for each value of $k$, and the second pass is only performed so far as to ensure a given fraction of the final optimal objective value. Large savings may be made in the second pass by choosing the required fraction to be significantly smaller than one.

A cheaper alternative is to use the Steihuag-Toint strategy, which is simply to stop at the first boundary point encountered along the piecewise linear path generated by the conjugate-gradient iterates. Note that if $H$ is significantly indefinite, this strategy often produces a far from optimal point, but is effective when $H$ is positive definite or almost so.

reference#

The method is described in detail in

N. I. M. Gould, S. Lucidi, M. Roma and Ph. L. Toint, ``Solving the trust-region subproblem using the Lanczos method’’. SIAM Journal on Optimization 9(2) (1999), 504-525.

introduction to function calls#

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

gltr_initialize - provide default control parameters and set up initial data structures
gltr_read_specfile (optional) - override control values by reading replacement values from a file
gltr_import_control - import control parameters prior to solution
gltr_solve_problem - solve the problem by reverse communication, a sequence of calls are made under control of a status parameter, each exit either asks the user to provide additional informaton and to re-enter, or reports that either the solution has been found or that an error has occurred
gltr_information (optional) - recover information about the solution and solution process
gltr_terminate - deallocate data structures

See the examples section for illustrations of use.

callable functions#

overview of functions provided#

// typedefs

typedef float spc_;
typedef double rpc_;
typedef int ipc_;

// structs

struct gltr_control_type;
struct gltr_inform_type;

// global functions

void gltr_initialize(
    void **data,
    struct gltr_control_type* control,
    ipc_ *status
);

void gltr_read_specfile(
    struct gltr_control_type* control,
    const char specfile[]
);

void gltr_import_control(
    struct gltr_control_type* control,
    void **data,
    ipc_ *status
);

void gltr_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ radius,
    rpc_ x[],
    rpc_ r[],
    rpc_ vector[]
);

void gltr_information(void **data, struct gltr_inform_type* inform, ipc_ *status);

void gltr_terminate(
    void **data,
    struct gltr_control_type* control,
    struct gltr_inform_type* inform
);

typedefs#

typedef float spc_

spc_ is real single precision

typedef double rpc_

rpc_ is the real working precision used, but may be changed to float by defining the preprocessor variable SINGLE.

typedef int ipc_

ipc_ is the default integer word length used, but may be changed to int64_t by defining the preprocessor variable INTEGER_64.

function calls#

void gltr_initialize(
    void **data,
    struct gltr_control_type* control,
    ipc_ *status
)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a struct containing control information (see gltr_control_type)

status

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

0
The initialization was successful.

void gltr_read_specfile(
    struct gltr_control_type* control,
    const char 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/gltr/GLTR.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/gltr.pdf for a list of how these keywords relate to the components of the control structure.

Parameters:

control	is a struct containing control information (see gltr_control_type)
specfile	is a character string containing the name of the specification file

void gltr_import_control(
    struct gltr_control_type* control,
    void **data,
    ipc_ *status
)

Import control parameters prior to solution.

Parameters:

control

is a struct whose members provide control paramters for the remaining prcedures (see gltr_control_type)

data

holds private internal data

status

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

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

void gltr_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ radius,
    rpc_ x[],
    rpc_ r[],
    rpc_ vector[]
)

Solve the trust-region problem using reverse communication.

Parameters:

data	holds private internal data
status	is a scalar variable of type ipc_, that gives the entry and exit status from the package. This must be set to 1 on initial entry. Set r (below) to $c$ for this entry. 4 the iteration is to be restarted with a smaller radius but with all other data unchanged. Set r (below) to $c$ for this entry. Possible exit values are: 0 the solution has been found 2 the inverse of $M$ must be applied to vector with the result returned in vector and the function re-entered with all other data unchanged. This will only happen if control.unitm is false 3 the product $H$ \* vector must be formed, with the result returned in vector and the function re-entered with all other data unchanged 5 The iteration must be restarted. Reset r (below) to $c$ and re-enter with all other data unchanged. This exit will only occur if control.steihaug_toint is false and the solution lies on the trust-region boundary -1 an array allocation has failed -2 an array deallocation has failed -3 n and/or radius is not positive -15 the matrix $M$ appears to be indefinite -18 the iteration limit has been exceeded -30 the trust-region has been encountered in Steihaug-Toint mode -31 the function value is smaller than control.f_min
n	is a scalar variable of type ipc_, that holds the number of variables
radius	is a scalar of type rpc_, that holds the trust-region radius, $\Delta$, used. radius must be strictly positive
x	is a one-dimensional array of size n and type rpc_, that holds the solution $x$. The j-th component of x, j = 0, … , n-1, contains $x_j$.
r	is a one-dimensional array of size n and type rpc_, that that must be set to $c$ on entry (status = 1) and re-entry ! (status = 4, 5). On exit, r contains the resiual $H x + c$.
vector	is a one-dimensional array of size n and type rpc_, that should be used and reset appropriately when status = 2 and 3 as directed.

void gltr_information(void **data, struct gltr_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see gltr_inform_type)

status

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

0
The values were recorded successfully

void gltr_terminate(
    void **data,
    struct gltr_control_type* control,
    struct gltr_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see gltr_control_type)
inform	is a struct containing output information (see gltr_inform_type)

available structures#

gltr_control_type structure#

#include <galahad_gltr.h>

struct gltr_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ itmax;
    ipc_ Lanczos_itmax;
    ipc_ extra_vectors;
    ipc_ ritz_printout_device;
    rpc_ stop_relative;
    rpc_ stop_absolute;
    rpc_ fraction_opt;
    rpc_ f_min;
    rpc_ rminvr_zero;
    rpc_ f_0;
    bool unitm;
    bool steihaug_toint;
    bool boundary;
    bool equality_problem;
    bool space_critical;
    bool deallocate_error_fatal;
    bool print_ritz_values;
    char ritz_file_name[31];
    char prefix[31];
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ error

error and warning diagnostics occur on stream error

ipc_ out

general output occurs on stream out

ipc_ print_level

the level of output required is specified by print_level

ipc_ itmax

the maximum number of iterations allowed (-ve = no bound)

ipc_ Lanczos_itmax

the maximum number of iterations allowed once the boundary has been encountered (-ve = no bound)

ipc_ extra_vectors

the number of extra work vectors of length n used

ipc_ ritz_printout_device

the unit number for writing debug Ritz values

rpc_ stop_relative

the iteration stops successfully when the gradient in the M(inverse) nor is smaller than max( stop_relative * initial M(inverse) gradient norm, stop_absolute )

rpc_ stop_absolute

see stop_relative

rpc_ fraction_opt

an estimate of the solution that gives at least .fraction_opt times the optimal objective value will be found

rpc_ f_min

the iteration stops if the objective-function value is lower than f_min

rpc_ rminvr_zero

the smallest value that the square of the M norm of the gradient of the the objective may be before it is considered to be zero

rpc_ f_0

the constant term, $f_0$, in the objective function

bool unitm

is $M$ the identity matrix ?

bool steihaug_toint

should the iteration stop when the Trust-region is first encountered ?

bool boundary

is the solution thought to lie on the constraint boundary ?

bool equality_problem

is the solution required to lie on the constraint boundary ?

bool space_critical

if .space_critical true, every effort will be made to use as little space as possible. This may result in longer computation time

bool deallocate_error_fatal

if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue

bool print_ritz_values

should the Ritz values be written to the debug stream?

char ritz_file_name[31]

name of debug file containing the Ritz values

char prefix[31]

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’

gltr_inform_type structure#

#include <galahad_gltr.h>

struct gltr_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ iter;
    ipc_ iter_pass2;
    rpc_ obj;
    rpc_ multiplier;
    rpc_ mnormx;
    rpc_ piv;
    rpc_ curv;
    rpc_ rayleigh;
    rpc_ leftmost;
    bool negative_curvature;
    bool hard_case;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

return status. See gltr_solve_problem for details

ipc_ alloc_status

the status of the last attempted allocation/deallocation

char bad_alloc[81]

the name of the array for which an allocation/deallocation error occurred

ipc_ iter

the total number of iterations required

ipc_ iter_pass2

the total number of pass-2 iterations required if the solution lies on the trust-region boundary

rpc_ obj

the value of the quadratic function

rpc_ multiplier

the Lagrange multiplier corresponding to the trust-region constraint

rpc_ mnormx

the $M$ -norm of $x$

rpc_ piv

the latest pivot in the Cholesky factorization of the Lanczos tridiagona

rpc_ curv

the most negative cuurvature encountered

rpc_ rayleigh

the current Rayleigh quotient

rpc_ leftmost

an estimate of the leftmost generalized eigenvalue of the pencil $(H,M)$

bool negative_curvature

was negative curvature encountered ?

bool hard_case

did the hard case occur ?

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/gltr/C/gltrt.c .

The floating-point type rpc_ is set in galahad_precision.h to double by default, but to float if the preprocessor variable SINGLE is defined. Similarly, the integer type ipc_ from galahad_precision.h is set to int by default, but to int64_t if the preprocessor variable INTEGER_64 is defined.

/* gltrt.c */
/* Full test for the GLTR C interface */

#include <stdio.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_gltr.h"

int main(void) {

    // Derived types
    void *data;
    struct gltr_control_type control;
    struct gltr_inform_type inform;

    // Set problem data
    ipc_ n = 100; // dimension

    ipc_ status;
    rpc_ radius;
    rpc_ x[n];
    rpc_ r[n];
    rpc_ vector[n];
    rpc_ h_vector[n];

    // Initialize gltr
    gltr_initialize( &data, &control, &status );

    // use a unit M ?
    for( ipc_ unit_m=0; unit_m <= 1; unit_m++){
      if ( unit_m == 0 ){
        control.unitm = false;
      } else {
        control.unitm = true;
      }
      gltr_import_control( &control, &data, &status );

      // resolve with a smaller radius ?
      for( ipc_ new_radius=0; new_radius <= 1; new_radius++){
        if ( new_radius == 0 ){
           radius = 1.0;
           status = 1;
        } else {
           radius = 0.1;
           status = 4;
        }
        for( ipc_ i = 0; i < n; i++) r[i] = 1.0;

        // iteration loop to find the minimizer
        while(true){ // reverse-communication loop
          gltr_solve_problem( &data, &status, n, radius, x, r, vector );
          if ( status == 0 ) { // successful termination
              break;
          } else if ( status < 0 ) { // error exit
              break;
          } else if ( status == 2 ) { // form the preconditioned vector
            for( ipc_ i = 0; i < n; i++) vector[i] = vector[i] / 2.0;
          } else if ( status == 3 ) { // form the Hessian-vector product
            h_vector[0] =  2.0 * vector[0] + vector[1];
            for( ipc_ i = 1; i < n-1; i++){
              h_vector[i] = vector[i-1] + 2.0 * vector[i] + vector[i+1];
            }
            h_vector[n-1] = vector[n-2] + 2.0 * vector[n-1];
            for( ipc_ i = 0; i < n; i++) vector[i] = h_vector[i];
          } else if ( status == 5 ) { // restart
            for( ipc_ i = 0; i < n; i++) r[i] = 1.0;
          }else{
              printf(" the value %1" i_ipc_ " of status should not occur\n",
                status);
              break;
          }
        }
        gltr_information( &data, &inform, &status );
        printf("MR = %1" i_ipc_ "%1" i_ipc_
               " gltr_solve_problem exit status = %" i_ipc_ ", f = %.2f\n",
               unit_m, new_radius, inform.status, inform.obj );
      }
    }
   // Delete internal workspace
   gltr_terminate( &data, &control, &inform );
}