GALAHAD LSRT package#

purpose#

The lsrt package uses a Krylov-subspace iteration to find an approximation of the global minimizer of the regularized linear sum-of-squares objective function. The aim is to minimize the least-squares objective function

\[r(x) = \frac{1}{2}\|Ax - b\|_2^2 + \frac{\sigma}{p} \|x\|_2^p,\]

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

method#

The required solution $x$ necessarily satisfies the optimality condition $A^T ( A x - b ) + \lambda x = 0$, where $\lambda = \sigma \|x\|_2^{p-2}.$

\noindent 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\|_2 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 regularized least-squares problem

\[\min \frac{1}{2} \| B_k y - \|b\| e_1 \|_2^2 + \frac{1}{p} \sigma \| y \|_2^p.\;\;\mbox{(1)}\]

To minimize (1), the optimality conditions

\[( B_k^T ( B_k^{} y(\lambda) - \|b\| e_1^{} ) + \lambda y(\lambda) = 0,\]

where $\lambda = \sigma \|y(\lambda)\|_2^{p-2}$, are used as the basis of an iteration. 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}\]

Thus, given an estimate $\lambda \geq 0$, this regularized least-squares problem 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) \|_2^{p-2} - \frac{\lambda}{\sigma} = 0.\]

In practice this nonlinear equation is reformulated, and a more rapidly converging iteration is used. Having found $y_k$, 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.

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

references#

A complete description of the un- an quadratically-regularized cases 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 the Newton-like process needed to determine $\lambda$ and other details are described 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 lsrt package must be called in the following order:

lsrt_initialize - provide default control parameters and set up initial data structures
lsrt_read_specfile (optional) - override control values by reading replacement values from a file
lsrt_import_control - import control parameters prior to solution
lsrt_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
lsrt_information (optional) - recover information about the solution and solution process
lsrt_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 lsrt_control_type;
struct lsrt_inform_type;

// global functions

void lsrt_initialize(
    void **data,
    struct lsrt_control_type* control,
    ipc_ *status
);

void lsrt_read_specfile(
    struct lsrt_control_type* control,
    const char specfile[]
);

void lsrt_import_control(
    struct lsrt_control_type* control,
    void **data,
    ipc_ *status
);

void lsrt_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const rpc_ power,
    const rpc_ weight,
    rpc_ x[],
    rpc_ u[],
    rpc_ v[]
);

void lsrt_information(void **data, struct lsrt_inform_type* inform, ipc_ *status);

void lsrt_terminate(
    void **data,
    struct lsrt_control_type* control,
    struct lsrt_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 calls#

void lsrt_initialize(
    void **data,
    struct lsrt_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 lsrt_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 lsrt_read_specfile(
    struct lsrt_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/lsrt/LSRT.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/lsrt.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 lsrt_control_type)
specfile	is a character string containing the name of the specification file

void lsrt_import_control(
    struct lsrt_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 lsrt_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 lsrt_solve_problem(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const rpc_ power,
    const rpc_ weight,
    rpc_ x[],
    rpc_ u[],
    rpc_ v[]
)

Solve the regularized least-squuares 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. 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, power 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$
power	is a scalar of type rpc_, that holds the regularization power, $p \geq 2$
weight	is a scalar of type rpc_, that holds the regularization weight, $\sigma > 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 lsrt_information(void **data, struct lsrt_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see lsrt_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 lsrt_terminate(
    void **data,
    struct lsrt_control_type* control,
    struct lsrt_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see lsrt_control_type)
inform	is a struct containing output information (see lsrt_inform_type)

available structures#

lsrt_control_type structure#

#include <galahad_lsrt.h>

struct lsrt_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_ bitmax;
    ipc_ extra_vectors;
    ipc_ stopping_rule;
    ipc_ freq;
    rpc_ stop_relative;
    rpc_ stop_absolute;
    rpc_ fraction_opt;
    rpc_ time_limit;
    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_ bitmax

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

ipc_ extra_vectors

the number of extra work vectors of length n used

ipc_ stopping_rule

the stopping rule used: 0=1.0, 1=norm step, 2=norm step/sigma (NOT USED)

ipc_ freq

frequency for solving the reduced tri-diagonal problem (NOT 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 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’

lsrt_inform_type structure#

#include <galahad_lsrt.h>

struct lsrt_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_ obj;
    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 lsrt_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

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 largest number of inner iterations performed during an outer iteration

rpc_ obj

the value of the objective function

rpc_ multiplier

the multiplier, $\lambda = sigma ||x||^(p-2)$

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 regularized linear least-squares problem; the code is available in $GALAHAD/src/lsrt/C/lsrtt.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.

/* lsrtt.c */
/* Full test for the LSRT C interface */

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

int main(void) {

    // Derived types
    void *data;
    struct lsrt_control_type control;
    struct lsrt_inform_type inform;

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

    ipc_ status;
    rpc_ power = 3.0;
    rpc_ weight = 1.0;
    rpc_ x[n];
    rpc_ u[m];
    rpc_ v[n];

    // Initialize lsrt
    lsrt_initialize( &data, &control, &status );

    status = 1;
    control.print_level = 0;
    lsrt_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
      lsrt_solve_problem( &data, &status, m, n, power, weight, 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;
      }
    }
    lsrt_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_lsrt.h
#include "galahad_pquad_lsrt.h"
#else
    printf("lsrt_solve_problem exit status = %" i_ipc_
           ", f = %.2f\n", inform.status, inform.obj );
    // Delete internal workspace
    lsrt_terminate( &data, &control, &inform );
#endif
}