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 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 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"

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 );
      printf("%1" i_ipc_ " lstr_solve_problem exit status = %" i_ipc_
             ", f = %.2f\n", new_radius, inform.status, inform.r_norm );
    }
   // Delete internal workspace
   lstr_terminate( &data, &control, &inform );
}