GALAHAD L2RT package#

purpose#

The l2rt package uses a Krylov-subspace iteration to find an approximation of the global minimizer of the regularized linear Euclidean-norm objective function. The aim is to minimize the objective function

\[r(x) = \sqrt{\|Ax - b\|_2^2 + \mu \|x\|_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 shift $\mu \geq 0$, 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/l2rt.pdf for additional details.

method#

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

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 Euclidean-norm problem

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

To solve (1), the optimality conditions

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

where

\[\lambda = \mu + \sigma \|y(\lambda)\|_2^{p-2} \sqrt{ \| B_k y(\lambda) - \|b\| e_1 \|_2^2 + \mu \|y(\lambda)\|_2^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 \sigma \|y(\lambda)\|_2^{p-2} \sqrt{ \| B_k y(\lambda) - \|b\| e_1 \|_2^2 + \mu \|y(\lambda)\|_2^2} + \mu - \lambda = 0.\;\;\mbox{(2)}\]

In practice (2) 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 the equation (2) significantly simplifies.

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 l2rt package must be called in the following order:

l2rt_initialize - provide default control parameters and set up initial data structures
l2rt_read_specfile (optional) - override control values by reading replacement values from a file
l2rt_import_control - import control parameters prior to solution
l2rt_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
l2rt_information (optional) - recover information about the solution and solution process
l2rt_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 l2rt_control_type;
struct l2rt_inform_type;

// global functions

void l2rt_initialize(
    void **data,
    struct l2rt_control_type* control,
    ipc_ *status
);

void l2rt_read_specfile(
    struct l2rt_control_type* control,
    const char specfile[]
);

void l2rt_import_control(
    struct l2rt_control_type* control,
    void **data,
    ipc_ *status
);

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

void l2rt_information(void **data, struct l2rt_inform_type* inform, ipc_ *status);

void l2rt_terminate(
    void **data,
    struct l2rt_control_type* control,
    struct l2rt_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 l2rt_initialize below will instead be

void l2rt_initialize_s_64(void **data, struct l2rt_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 l2rt_initialize(
    void **data,
    struct l2rt_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 l2rt_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 l2rt_read_specfile(
    struct l2rt_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/l2rt/L2RT.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/l2rt.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 l2rt_control_type)
specfile	is a character string containing the name of the specification file

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

Solve the regularized-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. 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, weight or shift 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$
shift	is a scalar of type rpc_, that holds the shift, $\mu$
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 l2rt_information(void **data, struct l2rt_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see l2rt_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 l2rt_terminate(
    void **data,
    struct l2rt_control_type* control,
    struct l2rt_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see l2rt_control_type)
inform	is a struct containing output information (see l2rt_inform_type)

available structures#

l2rt_control_type structure#

#include <galahad_l2rt.h>

struct l2rt_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’

l2rt_inform_type structure#

#include <galahad_l2rt.h>

struct l2rt_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 l2rt_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 = \mu + \sigma \|x\|^{p-2} * \sqrt{\|Ax-b\|^2 + \mu \|x\|^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 least-Euclidean-norm subproblem; the code is available in $GALAHAD/src/l2rt/C/l2rtt.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.

/* l2rtt.c */
/* Full test for the L2RT C interface */

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

int main(void) {

    // Derived types
    void *data;
    struct l2rt_control_type control;
    struct l2rt_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_ shift = 1.0;
    rpc_ x[n];
    rpc_ u[m];
    rpc_ v[n];

    // Initialize l2rt
    l2rt_initialize( &data, &control, &status );

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