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.

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

    function lstr_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 lstr_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 lstr_solve_problem(T, data, status, m, n, radius, x, u, v)

Solve the trust-region least-squares 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. 5 the iteration is to be restarted with a smaller radius but with all other data unchanged. Set `u` 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 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$
radius	is a scalar of type T that holds the trust-region radius, $\Delta > 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 lstr_information(T, data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see lstr_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 lstr_terminate(T, data, control, inform)

Deallocate all internal private storage

Parameters:

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

available structures#

lstr_control_type structure#

    struct lstr_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
      itmax_on_boundary::Int32
      bitmax::Int32
      extra_vectors::Int32
      stop_relative::T
      stop_absolute::T
      fraction_opt::T
      time_limit::T
      steihaug_toint::Bool
      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 itmax_on_boundary

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

Int32 bitmax

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

Int32 extra_vectors

the number of extra work vectors of length n 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 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

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’

lstr_inform_type structure#

    struct lstr_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
      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 lstr_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 if the solution lies on the trust-region boundary

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

T multiplier

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

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 linear least-squares trust-region problem; the code is available in $GALAHAD/src/lstr/Julia/test_lstr.jl .

# test_lstr.jl
# Simple code to test the Julia interface to LSTR

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_lstr(::Type{T}) where T
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{lstr_control_type{T}}()
  inform = Ref{lstr_inform_type{T}}()

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

  status = Ref{Cint}()
  radius = Ref{T}()
  x = zeros(T, n)
  u = zeros(T, m)
  v = zeros(T, n)

  # Initialize lstr
  lstr_initialize(T, data, control, status)

  # resolve with a smaller radius ?
  for new_radius in 0:1
    if new_radius == 0 # original radius
      radius[] = 1.0
      status[] = 1
    else # smaller radius
      radius[] = 0.1
      status[] = 5
    end

    @reset control[].print_level = Cint(0)
    lstr_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
      lstr_solve_problem(T, data, status, m, n, radius[], 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

    lstr_information(T, data, inform, status)
    # @printf("%1i lstr_solve_problem exit status = %i, f = %.2f\n", new_radius,
    #         inform[].status, inform[].r_norm)
  end

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

  return 0
end

@testset "LSTR" begin
  @test test_lstr(Float32) == 0
  @test test_lstr(Float64) == 0
  @test test_lstr(Float128) == 0
end