GALAHAD GLRT package

purpose

The glrt package uses a Krylov-subspace iteration to find an approximation of the global minimizer of regularized quadratic objective function. The aim is to minimize the regularized quadratic objective function

\[r(x) = f + g^T x + \frac{1}{2} x^T H x + \frac{\sigma}{p} \|x\|_{M}^p,\]

where the weight \(\sigma \geq 0\), the power \(p \geq 2\), and where the \(M\)-norm of \(x\) is defined to be \(\|x\|_{M} = \sqrt{x^T M x}\). 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/glrt.pdf for additional details.

method

The required solution \(x\) necessarily satisfies the optimality condition \(H x + \lambda M x + g = 0\), where \(\lambda = \sigma \|x\|_{M}^{p-2}\). In addition, the matrix \(H + \lambda M\) will be positive semi-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 + \frac{1}{p} \sigma \| y \|_2^p,\;\;\mbox{(1)}\]

where \(e_1\) is the first unit vector.

To minimize (1), the optimality conditions

\[( T_k + \lambda I ) y(\lambda) = - g,\;\;\mbox{(2)}\]

where \(\lambda = \sigma \|y(\lambda)\|_{M}^{p-2}\) are used as the basis of an iteration. Specifically, given an estimate \(\lambda\) for which \(T_k + \lambda I\) is positive definite, the tridiagonal system (2) may be efficiently solved to give \(y(\lambda)\). It is then simply a matter of adjusting \(\lambda\) (for example by a Newton-like process) to solve the scalar nonlinear equation

\[\theta(\lambda) \equiv \|y(\lambda)\|_{M}^{p-2} - \frac{\lambda}{\sigma} = 0.\;\;\mbox{(3)}\]

In practice (3) is reformulated, and a more rapidly converging iteration is used.

It is possible to measure 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.

Special code is used in the special case \(p=2\), as in this case a single pass suffices.

reference

The method is described in detail in

C. Cartis, N. I. M. Gould and Ph. L. Toint, ``Adaptive cubic regularisation methods for unconstrained optimization. Part {I}: motivation, convergence and numerical results’’. Mathematical Programming 127(2) (2011) 245-295.

introduction to function calls

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

• glrt_initialize - provide default control parameters and set up initial data structures

• glrt_read_specfile (optional) - override control values by reading replacement values from a file

• glrt_import_control - import control parameters prior to solution

• glrt_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

• glrt_information (optional) - recover information about the solution and solution process

• glrt_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 glrt_control_type;
struct glrt_inform_type;

// global functions

void glrt_initialize(
    void **data,
    struct glrt_control_type* control,
    ipc_ *status
);

void glrt_read_specfile(
    struct glrt_control_type* control,
    const char specfile[]
);

void glrt_import_control(
    struct glrt_control_type* control,
    void **data,
    ipc_ *status
);

void glrt_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ power,
    const rpc_ weight,
    rpc_ x[],
    rpc_ r[],
    rpc_ vector[]
);

void glrt_information(void **data, struct glrt_inform_type* inform, ipc_ *status);

void glrt_terminate(
    void **data,
    struct glrt_control_type* control,
    struct glrt_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 glrt_initialize(
    void **data,
    struct glrt_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 glrt_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 glrt_read_specfile(
    struct glrt_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/glrt/GLRT.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/glrt.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 glrt_control_type)
specfile	is a character string containing the name of the specification file

void glrt_import_control(
    struct glrt_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 glrt_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 glrt_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ power,
    const rpc_ weight,
    rpc_ x[],
    rpc_ r[],
    rpc_ vector[]
)

Solve the regularized-quadratic 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. 6 the iteration is to be restarted with a larger weight 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 4** The iteration must be restarted. Reset r (below) to \(c\) and re-enter with all other data unchanged. -1 an array allocation has failed -2 an array deallocation has failed -3 n and/or radius is not positive -7 the problem is unbounded from below. This can only happen if power = 2, and in this case the objective is unbounded along the arc x + t vector as t goes to infinity -15 the matrix \(M\) appears to be indefinite -18 the iteration limit has been exceeded
n	is a scalar variable of type ipc_, that holds the number of variables
power	is a scalar of type rpc_, that holds the egularization power, \(p \geq 2\)
weight	is a scalar of type rpc_, that holds the positive regularization weight, \(\sigma\)
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 glrt_information(void **data, struct glrt_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data	holds private internal data
inform	is a struct containing output information (see glrt_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 glrt_terminate(
    void **data,
    struct glrt_control_type* control,
    struct glrt_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see glrt_control_type)
inform	is a struct containing output information (see glrt_inform_type)

available structures

glrt_control_type structure

#include <galahad_glrt.h>

struct glrt_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ itmax;
    ipc_ stopping_rule;
    ipc_ freq;
    ipc_ extra_vectors;
    ipc_ ritz_printout_device;
    rpc_ stop_relative;
    rpc_ stop_absolute;
    rpc_ fraction_opt;
    rpc_ rminvr_zero;
    rpc_ f_0;
    bool unitm;
    bool impose_descent;
    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_ stopping_rule

the stopping rule used (see below). Possible values are:

• 1 stopping rule = norm of the step.

• 2 stopping rule is norm of the step / \(\sigma\).

• other. stopping rule = 1.0.

ipc_ freq

frequency for solving the reduced tri-diagonal problem

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^{-1}\) norm is smaller than max( stop_relative \* min( 1, stopping_rule ) \* norm initial gradient, 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_ rminvr_zero

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

rpc_ f_0

the constant term, f0, in the objective function

bool unitm

is M the identity matrix ?

bool impose_descent

is descent required i.e., should \(c^T x < 0\)?

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’

glrt_inform_type structure

#include <galahad_glrt.h>

struct glrt_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ iter;
    ipc_ iter_pass2;
    rpc_ obj;
    rpc_ obj_regularized;
    rpc_ multiplier;
    rpc_ xpo_norm;
    rpc_ leftmost;
    bool negative_curvature;
    bool hard_case;
};

detailed documentation

inform derived type as a C struct

components

ipc_ status

return status. See glrt_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

rpc_ obj

the value of the quadratic function

rpc_ obj_regularized

the value of the regularized quadratic function

rpc_ multiplier

the multiplier, \(\sigma \|x\|^{p-2}\)

rpc_ xpo_norm

the value of the norm \(\|x\|_M\)

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 regularization subproblem; the code is available in $GALAHAD/src/glrt/C/glrtt.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.

/* glrtt.c */
/* Full test for the GLRT C interface */

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

int main(void) {

    // Derived types
    void *data;
    struct glrt_control_type control;
    struct glrt_inform_type inform;

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

    ipc_ status;
    rpc_ weight;
    rpc_ power = 3.0;
    rpc_ x[n];
    rpc_ r[n];
    rpc_ vector[n];
    rpc_ h_vector[n];

    // Initialize glrt
    glrt_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;
      }
      glrt_import_control( &control, &data, &status );
      // resolve with a larger weight ?
      for( ipc_ new_weight=0; new_weight <= 1; new_weight++){
        if ( new_weight == 0 ){
           weight = 1.0;
           status = 1;
        } else {
           weight = 10.0;
           status = 6;
        }
        for( ipc_ i = 0; i < n; i++) r[i] = 1.0;

        // iteration loop to find the minimizer
        while(true){ // reverse-communication loop
          glrt_solve_problem( &data, &status, n, power, weight, 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 == 4 ) { // 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;
          }
        }
        glrt_information( &data, &inform, &status );
        printf("MR = %1" i_ipc_ "%1" i_ipc_
               " glrt_solve_problem exit status = %" i_ipc_ ", f = %.2f\n",
               unit_m, new_weight, inform.status, inform.obj_regularized );
      }
    }
   // Delete internal workspace
   glrt_terminate( &data, &control, &inform );
}