GALAHAD LSTR package

purpose

The lstr package uses a Krylov-subspace iteration to find an approximation of the global minimizer of the linear sum-of-squares objective function within a sphere; this is commonly known as the linear least-squares trust-region subproblem. The aim is to minimize the least-squares objective function

\[q(x) = \frac{1}{2}\|Ax - b\|_2^2,\]

where the vector \(x\) is required to satisfy the spherical trust-region constraint \(\|x\|_2 \leq \Delta\), where the \(\ell_2\)-norm of \(x\) is defined to be \(\|x\|_2 = \sqrt{x^T x}\), and where the radius \(\Delta > 0\). The method may be suitable for large problems as no factorization is required. Reverse communication is used to obtain matrix-vector products of the form \(u + A v\) and \(v + A^T u.\)

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

method

The required solution \(x\) necessarily satisfies the optimality condition \(A^T ( A x - b ) + \lambda x = 0\), where \(\lambda \geq 0\) is a Lagrange multiplier corresponding to the trust-region constraint \(\|x\|_2 \leq \Delta\).

The method is iterative. Starting with the vector \(u_1 = b\), a bi-diagonalisation process is used to generate the vectors \(v_k\) and \(u_k+1\) so that the \(n\) by \(k\) matrix \(V_k = ( v_1 \ldots v_k)\) and the \(m\) by \((k+1)\) matrix \(U_k = ( u_1 \ldots u_{k+1})\) together satisfy

\[A V_k = U_{k+1} B_k \;\;\mbox{and}\;\; b = \|b\| U_{k+1} e_1,\]

where \(B_k\) is \((k+1)\) by \(k\) and lower bi-diagonal, \(U_k\) and \(V_k\) have orthonormal columns and \(e_1\) is the first unit vector. The solution sought is of the form \(x_k = V_k y_k\), where \(y_k\) solves the bi-diagonal least-squares trust-region problem

\[\min \| B_k y - \|b\| e_1 \|_2 \;\;\mbox{subject to}\;\; \|y\|_2 \leq \Delta.\;\;\mbox{(1)}\]

If the trust-region constraint is inactive, the solution \(y_k\) may be found, albeit indirectly, via the LSQR algorithm of Paige and Saunders which solves the bi-diagonal least-squares problem

\[\min \| B_k y - \|b\| e_1 \|_2\]

using a QR factorization of \(B_k\). Only the most recent \(v_k\) and \(u_{k+1}\) are required, and their predecessors discarded, to compute \(x_k\) from \(x_{k-1}\). This method has the important property that the iterates \(y\) (and thus \(x_k\)) generated increase in norm with \(k\). Thus as soon as an LSQR iterate lies outside the trust-region, the required solution to (1) and thus to the original problem must lie on the boundary of the trust-region.

If the solution is so constrained, the simplest strategy is to interpolate the last interior iterate with the newly discovered exterior one to find the boundary point — the so-called Steihaug-Toint point — between them. Once the solution is known to lie on the trust-region boundary, further improvement may be made by solving

\[\min \| B_k y - \|b\| e_1 \|_2 \;\;\mbox{subject to}\;\; \|y\|_2 = \Delta,\]

for which the optimality conditions require that \(y_k = y(\lambda_k)\) where \(\lambda_k\) is the positive root of

\[B_k^T ( B_k^{} y(\lambda) - \|b\| e_1^{} ) + \lambda y(\lambda) = 0 \;\;\mbox{and}\;\; \|y(\lambda)\|_2 = \Delta\]

The vector \(y(\lambda)\) is equivalently the solution to the regularized least-squares problem

\[\begin{split}\min \left \| \left(\begin{array}{c} B_k \\ \lambda^{1/2} I \end{array}\right) y - \|b\| e_1^{} \right \|_2\end{split}\]

and may be found efficiently. Given \(y(\lambda)\), Newton’s method is then used to find \(\lambda_k\) as the positive root of \(\|y(\lambda)\|_2 = \Delta\). Unfortunately, unlike when the solution lies in the interior of the trust-region, it is not known how to recur \(x_k\) from \(x_{k-1}\) given \(y_k\), and a second pass in which \(x_k = V_k y_k\) is regenerated is needed — this need only be done once \(x_k\) has implicitly deemed to be sufficiently close to optimality. 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.

references

A complete description of the unconstrained case is given by

C. C. Paige and M. A. Saunders, ``LSQR: an algorithm for sparse linear equations and sparse least squares’’. ACM Transactions on Mathematical Software 8(1 (1982) 43–71,

and

C. C. Paige and M. A. Saunders, ``ALGORITHM 583: LSQR: an algorithm for sparse linear equations and sparse least squares’’. ACM Transactions on Mathematical Software 8(2) (1982) 195–209.

Additional details on how to proceed once the trust-region constraint is encountered are described in detail in

C. Cartis, N. I. M. Gould and Ph. L. Toint, ``Trust-region and other regularisation of linear least-squares problems’’. BIT 49(1) (2009) 21-53.

introduction to function calls

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

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

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

• lstr_import_control - import control parameters prior to solution

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

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

• lstr_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 lstr_control_type;
struct lstr_inform_type;

// global functions

void lstr_initialize(
    void **data,
    struct lstr_control_type* control,
    ipc_ *status
);

void lstr_read_specfile(
    struct lstr_control_type* control,
    const char specfile[]
);

void lstr_import_control(
    struct lstr_control_type* control,
    void **data,
    ipc_ *status
);

void lstr_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const rpc_ radius,
    rpc_ x[],
    rpc_ u[],
    rpc_ v[]
);

void lstr_information(void **data, struct lstr_inform_type* inform, ipc_ *status);

void lstr_terminate(
    void **data,
    struct lstr_control_type* control,
    struct lstr_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 REAL_32 or (if supported) to __real128 using the variable REAL_128.

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 and structure names

The function and structure names described below are appropriate for the default real working precision (double) and integer word length (int32_t). To use the functions and structures with different precisions and integer word lengths, an additional suffix must be added to their names (and the arguments set accordingly). The appropriate suffices are:

_s for single precision (float) reals and standard 32-bit (int32_t) integers;

_q for quadruple precision (__real128) reals (if supported) and standard 32-bit (int32_t) integers;

_64 for standard precision (double) reals and 64-bit (int64_t) integers;

_s_64 for single precision (float) reals and 64-bit (int64_t) integers; and

_q_64 for quadruple precision (__real128) reals (if supported) and 64-bit (int64_t) integers.

Thus a call to lstr_initialize below will instead be

void lstr_initialize_s_64(void **data, struct lstr_control_type_s_64* control,
                         int64_t *status)

if single precision (float) reals and 64-bit (int64_t) integers are required. Thus it is possible to call functions for this package with more that one precision and/or integer word length at same time. An example is provided for the package expo, and the obvious modifications apply equally here.

function calls

void lstr_initialize(
    void **data,
    struct lstr_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 lstr_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 lstr_read_specfile(
    struct lstr_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/lstr/LSTR.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/lstr.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 lstr_control_type)
specfile	is a character string containing the name of the specification file

void lstr_import_control(
    struct lstr_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 lstr_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 lstr_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const rpc_ radius,
    rpc_ x[],
    rpc_ u[],
    rpc_ v[]
)

Solve the trust-region least-squares 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 u (below) to \(b\) for this entry. 5 the iteration is to be restarted with a smaller radius but with all other data unchanged. Set u (below) to \(b\) for this entry. Possible exit values are: 0 the solution has been found 2 The user must perform the operation \[u := u + A v,\] n u := u + A v, \n and recall the function. The vectors \(u\) and \(v\) are available in the arrays u and v (below) respectively, and the result \(u\) must overwrite the content of u. No argument except u should be altered before recalling the function 3 The user must perform the operation \[v := v + A^T u,\] n v := v + A^T u, \n and recall the function. The vectors \(u\) and \(v\) are available in the arrays u and v (below) respectively, and the result \(v\) must overwrite the content of v. No argument except v should be altered before recalling the function 4 The user must reset u (below) to \(b\) are recall the function. No argument except u should be altered before recalling the function -1 an array allocation has failed -2 an array deallocation has failed -3 one or more of n, m or weight violates allowed bounds -18 the iteration limit has been exceeded -25 status is negative on entry
m	is a scalar variable of type ipc_, that holds the number of equations (i.e., rows of \(A\)), \(m > 0\)
n	is a scalar variable of type ipc_, that holds the number of variables (i.e., columns of \(A\)), \(n > 0\)
radius	is a scalar of type rpc_, that holds the trust-region radius, \(\Delta > 0\)
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\).
u	is a one-dimensional array of size m and type rpc_, that should be used and reset appropriately when status = 1 to 5 as directed by status.
v	is a one-dimensional array of size n and type rpc_, that should be used and reset appropriately when status = 1 to 5 as directed by status.

void lstr_information(void **data, struct lstr_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data	holds private internal data
inform	is a struct containing output information (see lstr_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 lstr_terminate(
    void **data,
    struct lstr_control_type* control,
    struct lstr_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see lstr_control_type)
inform	is a struct containing output information (see lstr_inform_type)

available structures

lstr_control_type structure

#include <galahad_lstr.h>

struct lstr_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ start_print;
    ipc_ stop_print;
    ipc_ print_gap;
    ipc_ itmin;
    ipc_ itmax;
    ipc_ itmax_on_boundary;
    ipc_ bitmax;
    ipc_ extra_vectors;
    rpc_ stop_relative;
    rpc_ stop_absolute;
    rpc_ fraction_opt;
    rpc_ time_limit;
    bool steihaug_toint;
    bool space_critical;
    bool deallocate_error_fatal;
    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_ start_print

any printing will start on this iteration

ipc_ stop_print

any printing will stop on this iteration

ipc_ print_gap

the number of iterations between printing

ipc_ itmin

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

ipc_ itmax

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

ipc_ itmax_on_boundary

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

ipc_ bitmax

the maximum number of Newton inner iterations per outer iteration allowe (-ve = no bound)

ipc_ extra_vectors

the number of extra work vectors of length n used

rpc_ stop_relative

the iteration stops successfully when \(\|A^Tr\|\) is less than max( stop_relative \* \(\|A^Tr_{initial} \|\), 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_ time_limit

the maximum elapsed time allowed (-ve means infinite)

bool steihaug_toint

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

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

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’

lstr_inform_type structure

#include <galahad_lstr.h>

struct lstr_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ iter;
    ipc_ iter_pass2;
    ipc_ biters;
    ipc_ biter_min;
    ipc_ biter_max;
    rpc_ multiplier;
    rpc_ x_norm;
    rpc_ r_norm;
    rpc_ Atr_norm;
    rpc_ biter_mean;
};

detailed documentation

inform derived type as a C struct

components

ipc_ status

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

ipc_ biters

the total number of inner iterations performed

ipc_ biter_min

the smallest number of inner iterations performed during an outer iteration

ipc_ biter_max

the largestt number of inner iterations performed during an outer iteration

rpc_ multiplier

the Lagrange multiplier, \(\lambda\), corresponding to the trust-region constraint

rpc_ x_norm

the Euclidean norm of \(x\)

rpc_ r_norm

the Euclidean norm of \(Ax-b\)

rpc_ Atr_norm

the Euclidean norm of \(A^T (Ax-b) + \lambda x\)

rpc_ biter_mean

the average number of inner iterations performed during an outer

iteration

Index

example calls

This is an example of how to use the package to solve a linear least-squares trust-region problem; the code is available in $GALAHAD/src/lstr/C/lstrt.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.

/* lstrt.c */
/* Full test for the LSTR C interface */

#include <stdio.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_lstr.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

int main(void) {

    // Derived types
    void *data;
    struct lstr_control_type control;
    struct lstr_inform_type inform;

    // Set problem data
    ipc_ n = 50; // dimensions
    ipc_ m = 2 * n;

    ipc_ status;
    rpc_ radius;
    rpc_ x[n];
    rpc_ u[m];
    rpc_ v[n];

    // Initialize lstr
    lstr_initialize( &data, &control, &status );

    // resolve with a smaller radius ?
    for( ipc_ new_radius=0; new_radius <= 1; new_radius++){
      if ( new_radius == 0 ){ // original radius
         radius = 1.0;
         status = 1;
      } else { // smaller radius
         radius = 0.1;
         status = 5;
      }
      control.print_level = 0;
      lstr_import_control( &control, &data, &status );

      for( ipc_ i = 0; i < m; i++) u[i] = 1.0; // b = 1

      // iteration loop to find the minimizer with A^T = (I:diag(1:n))
      while(true){ // reverse-communication loop
        lstr_solve_problem( &data, &status, m, n, radius, x, u, v );
        if ( status == 0 ) { // successful termination
            break;
        } else if ( status < 0 ) { // error exit
            break;
        } else if ( status == 2 ) { // form u <- u + A * v
          for( ipc_ i = 0; i < n; i++) {
            u[i] = u[i] + v[i];
            u[n+i] = u[n+i] + (i+1)*v[i];
          }
        } else if ( status == 3 ) { // form v <- v + A^T * u
          for( ipc_ i = 0; i < n; i++) v[i] = v[i] + u[i] + (i+1) * u[n+i];
        } else if ( status == 4 ) { // restart
          for( ipc_ i = 0; i < m; i++) u[i] = 1.0;
        }else{
            printf(" the value %1" i_ipc_ " of status should not occur\n",
              status);
            break;
        }
      }
      lstr_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_lstr.h
#include "galahad_pquad_lstr.h"
#else
      printf("%1" i_ipc_ " lstr_solve_problem exit status = %" i_ipc_
             ", f = %.2f\n", new_radius, inform.status, inform.r_norm );
#endif
    }
   // Delete internal workspace
   lstr_terminate( &data, &control, &inform );
}