LSTR#

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.

functions#

lstr.initialize()#

Set default option values and initialize private data

Returns:

optionsdict

dictionary containing default control options:

errorint
error and warning diagnostics occur on stream error.

outint
general output occurs on stream out.

problemint
unit to write problem data into file problem_file.

print_levelint
the level of output required is specified by print_level. Possible values are

<=0

gives no output,

1

gives a one-line summary for every iteration.

2

gives a summary of the inner iteration for each iteration.

>=3

gives increasingly verbose (debugging) output.

start_printint
any printing will start on this iteration.

stop_printint
any printing will stop on this iteration.

print_gapint
the number of iterations between printing.

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

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

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

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

extra_vectorsint
the number of extra work vectors of length n used.

stop_relativefloat
the iteration stops successfully when $\|A^Tr\|$ is less than max( stop_relative * $\|A^T b \|$, stop_absolute ), where $r = A x - b$.

stop_absolutefloat
see stop_relative.

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

time_limitfloat
the maximum elapsed time allowed (-ve means infinite).

steihaug_tointbool
should the iteration stop when the Trust-region is first encountered?.

space_criticalbool
if space_critical is True, every effort will be made to use as little space as possible. This may result in longer computation time.

deallocate_error_fatalbool
if deallocate_error_fatal is True, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue.

prefixstr
all output lines will be prefixed by the string contained in quotes within prefix, e.g. ‘word’ (note the qutoes) will result in the prefix word.

lstr.load_options(options=None)#

Import control options into internal storage prior to solution.

Parameters:

optionsdict, optional
dictionary of control options (see lstr.initialize).

lstr.solve_problem(status, m, n, radius, u, v)#

Find the global moinimizer of the quadratic objective function $q(x)$ within a trust-region of radius $\Delta$.

Parameters:

statusint
holds the entry status. Possible values are

1

an initial entry with u set to $b$.

5

a restart entry with u reset to $b$, but a smaller radius $\Delta$.

other

the value returned from the previous call, see Returns below.

mint
holds the number of residuals, i.e., the number of rows of $A$.

nint
holds the number of variables, i.e., the number of columns of $A$.

radiusfloat
holds the strictly positive trust-region radius, $\Delta$.

undarray(m)
holds the result vector when initial or return status = 1, 2, 4 or 5 (see below).

vndarray(n)
holds the result vector when return status = 3 (see below).

Returns:

statusint
holds the exit status. Possible values are

0

the solution has been found, no further reentry is required

2

the sum $u + A v$, involving the vectors $u$ and $v$ returned in u and v, must be formed, the result placed in u, and the function recalled with status set to 2.

3

the sum $v + A^T u$, involving the vectors $u$ and $v$ returned in u and v, must be formed, the result placed in v, and the function recalled with status set to 3.

4

the iteration must be restarted by setting u to $b$, and the function recalled with status set to 4.

<0

an error occurred, see status in lstr.information for further details.

xndarray(n)
holds the values of the approximate minimizer $x$.

undarray(m)
holds the result vector $u$.

vndarray(n)
holds the result vector $v$.

[optional] lstr.information()

Provide optional output information

Returns:

informdict

dictionary containing output information:

statusint
return status. Possible values are:

0

The run was successful.

-1

An allocation error occurred. A message indicating the offending array is written on unit options[‘error’], and the returned allocation status and a string containing the name of the offending array are held in inform[‘alloc_status’] and inform[‘bad_alloc’] respectively.

-2

A deallocation error occurred. A message indicating the offending array is written on unit options[‘error’] and the returned allocation status and a string containing the name of the offending array are held in inform[‘alloc_status’] and inform[‘bad_alloc’] respectively.

-3

The restriction n > 0, m > 0 or $\Delta > 0$ has been violated.

-18

The iteration limit has been exceeded.

-25

status is negative on entry.

alloc_statusint
the status of the last attempted allocation/deallocation.

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

iterint
the total number of iterations required.

iter_pass2int
the total number of pass-2 iterations required if the solution lies on the trust-region boundary.

bitersint
the total number of inner iterations performed.

biter_minint
the smallest number of inner iterations performed during an outer iteration.

biter_maxint
the largestt number of inner iterations performed during an outer iteration.

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

x_normfloat
the Euclidean norm of $x$.

r_normfloat
the Euclidean norm of $Ax-b$.

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

biter_meanfloat
the average number of inner iterations performed during an outer iteration.

lstr.terminate()#

Deallocate all internal private storage.

example code#

from galahad import lstr
import numpy as np
np.set_printoptions(precision=4,suppress=True,floatmode='fixed')
print("\n** python test: lstr")

# set problem data

n = 50
m = 2 * n
u = np.empty(m)
v = np.empty(n)
b = np.empty(m)
b[:] = 1.0

# set default options
options = lstr.initialize()

# set some non-default options
options['unitm'] = False
options['print_level'] = 0
#print("options:", options)

# load non-default options
lstr.load_options(options)

# first problem

print(" solve problem 1")
status = 1
radius = 1.0
u = b
while status > 0:
  status, x, u, v = lstr.solve_problem(status, m, n, radius, u, v)
  if status == 2: # u -> u + A v
    for i in range(n):
      u[i] = u[i] + v[i]
      u[n+i] = u[n+i] + (i+1) * v[i]
  elif status == 3: # v -> v + A^T u
    for i in range(n):
      v[i] = v[i] + u[i] + (i+1) * u[n+i]
  elif status == 4: # reset r to b
    u = b

# get information

inform = lstr.information()
print(" optimal ||r||: %.4f" % inform['r_norm'])
#print(" x:",x)

# second problem, same data with smaller radius

print("\n solve problem 2 with smaller radius")
status = 5
radius = 0.1
u = b
while status > 0:
  status, x, u, v = lstr.solve_problem(status, m, n, radius, u, v)
  if status == 2: # u -> u + A v
    for i in range(n):
      u[i] = u[i] + v[i]
      u[n+i] = u[n+i] + (i+1) * v[i]
  elif status == 3: # v -> v + A^T u
    for i in range(n):
      v[i] = v[i] + u[i] + (i+1) * u[n+i]
  elif status == 4: # reset r to b
    u = b

# get information

inform = lstr.information()
print(" optimal ||r||: %.4f" % inform['r_norm'])
#print(" x:",x)
print('** lstr exit status:', inform['status'])

# deallocate internal data

lstr.terminate()

This example code is available in $GALAHAD/src/lstr/Python/test_lstr.py .