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.

parametric real type T#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision), Float64 (double precision) or, if supported, Float128 (quadruple precision).

callable functions#

    function lsrt_initialize(T, data, control, status)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a structure containing control information (see lsrt_control_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

0
The initialization was successful.

    function lsrt_read_specfile(T, control, 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 structure containing control information (see lsrt_control_type)
specfile	is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file

    function lsrt_import_control(T, control, data, status)

Import control parameters prior to solution.

Parameters:

control

is a structure whose members provide control parameters for the remaining procedures (see lsrt_control_type)

data

holds private internal data

status

is a scalar variable of type Int32 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

    function lsrt_solve_problem(T, data, status, m, n, power, weight, x, u, v)

Solve the regularized least-squuares problem using reverse communication.

Parameters:

data	holds private internal data
status	is a scalar variable of type Int32 that gives the entry and exit status from the package. This must be set to 1 on initial entry. Set the argument `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 + Av,\] 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^Tu,\] and recall the function. The vectors $u$ and $v$ are available in the arrays `u` and `v` 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` 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 Int32 that holds the number of equations (i.e., rows of $A$), $m > 0$
n	is a scalar variable of type Int32 that holds the number of variables (i.e., columns of $A$), $n > 0$
power	is a scalar of type T that holds the regularization power, $p \geq 2$
weight	is a scalar of type T that holds the regularization weight, $\sigma > 0$
x	is a one-dimensional array of size n and type T that holds the solution $x$. The j-th component of `x`, j = 1, … , n, contains $x_j$.
u	is a one-dimensional array of size m and type T 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 T that should be used and reset appropriately when status = 1 to 5 as directed by status.

    function lsrt_information(T, data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see lsrt_inform_type)

status

is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):

0
The values were recorded successfully

    function lsrt_terminate(T, data, control, inform)

Deallocate all internal private storage

Parameters:

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

available structures#

lsrt_control_type structure#

    struct lsrt_control_type{T}
      f_indexing::Bool
      error::Int32
      out::Int32
      print_level::Int32
      start_print::Int32
      stop_print::Int32
      print_gap::Int32
      itmin::Int32
      itmax::Int32
      bitmax::Int32
      extra_vectors::Int32
      stopping_rule::Int32
      freq::Int32
      stop_relative::T
      stop_absolute::T
      fraction_opt::T
      time_limit::T
      space_critical::Bool
      deallocate_error_fatal::Bool
      prefix::NTuple{31,Cchar}

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

Int32 error

error and warning diagnostics occur on stream error

Int32 out

general output occurs on stream out

Int32 print_level

the level of output required is specified by print_level

Int32 start_print

any printing will start on this iteration

Int32 stop_print

any printing will stop on this iteration

Int32 print_gap

the number of iterations between printing

Int32 itmin

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

Int32 itmax

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

Int32 bitmax

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

Int32 extra_vectors

the number of extra work vectors of length n used

Int32 stopping_rule

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

Int32 freq

frequency for solving the reduced tri-diagonal problem (NOT USED)

T stop_relative

the iteration stops successfully when ||A^Tr|| is less than max( stop_relative * ||A^Tr initial ||, stop_absolute )

T stop_absolute

see stop_relative

T fraction_opt

an estimate of the solution that gives at least .fraction_opt times the optimal objective value will be found

T 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

NTuple{31,Cchar} prefix

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#

    struct lsrt_inform_type{T}
      status::Int32
      alloc_status::Int32
      bad_alloc::NTuple{81,Cchar}
      iter::Int32
      iter_pass2::Int32
      biters::Int32
      biter_min::Int32
      biter_max::Int32
      obj::T
      multiplier::T
      x_norm::T
      r_norm::T
      Atr_norm::T
      biter_mean::T

detailed documentation#

inform derived type as a Julia structure

components#

Int32 status

return status. See lsrt_solve_problem for details

Int32 alloc_status

the status of the last attempted allocation/deallocation

NTuple{81,Cchar} bad_alloc

the name of the array for which an allocation/deallocation error occurred

Int32 iter

the total number of iterations required

Int32 iter_pass2

the total number of pass-2 iterations required

Int32 biters

the total number of inner iterations performed

Int32 biter_min

the smallest number of inner iterations performed during an outer iteration

Int32 biter_max

the largest number of inner iterations performed during an outer iteration

T obj

the value of the objective function

T multiplier

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

T x_norm

the Euclidean norm of $x$

T r_norm

the Euclidean norm of $Ax-b$

T Atr_norm

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

T 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/Julia/test_lsrt.jl .

# test_lsrt.jl
# Simple code to test the Julia interface to LSRT

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_lsrt(::Type{T}) where T
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{lsrt_control_type{T}}()
  inform = Ref{lsrt_inform_type{T}}()

  # Set problem data
  n = 50 # dimensions
  m = 2 * n

  status = Ref{Cint}()
  power = T(3.0)
  weight = one(T)
  x = zeros(T, n)
  u = zeros(T, m)
  v = zeros(T, n)

  # Initialize lsrt
  lsrt_initialize(T, data, control, status)

  status[] = 1
  @reset control[].print_level = Cint(0)
  lsrt_import_control(T, control, data, status)

  for i in 1:m
    u[i] = 1.0 # b = 1
  end

  # iteration loop to find the minimizer with A^T = (I:diag(1:n))
  terminated = false
  while !terminated # reverse-communication loop
    lsrt_solve_problem(T, data, status, m, n, power, weight, x, u, v)
    if status[] == 0 # successful termination
      terminated = true
    elseif status[] < 0 # error exit
      terminated = true
    elseif status[] == 2 # form u <- u + A * v
      for i in 1:n
        u[i] = u[i] + v[i]
        u[n + i] = u[n + i] + (i + 1) * v[i]
      end
    elseif status[] == 3 # form v <- v + A^T * u
      for i in 1:n
        v[i] = v[i] + u[i] + (i + 1) * u[n + i]
      end
    elseif status[] == 4 # restart
      for i in 1:m
        u[i] = 1.0
      end
    else
      @printf(" the value %1i of status should not occur\n",
              status)
    end
  end

  lsrt_information(T, data, inform, status)
  @printf("lsrt_solve_problem exit status = %i, f = %.2f\n", inform[].status, inform[].obj)

  # Delete internal workspace
  lsrt_terminate(T, data, control, inform)

  return 0
end

@testset "LSRT" begin
  @test test_lsrt(Float32) == 0
  @test test_lsrt(Float64) == 0
  @test test_lsrt(Float128) == 0
end