LSRT#

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.

functions#

lsrt.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).

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

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

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

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

stop_relativefloat
the iteration stops successfully when $\|A^Tr\|$ is less than max( stop_relative * $\|A^Tb\|$, 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).

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.

lsrt.load_options(options=None)#

Import control options into internal storage prior to solution.

Parameters:

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

lsrt.solve_problem(status, m, n, power, weight, u, v)#

Find the global moinimizer of the regularized quadratic objective function $r(x)$.

Parameters:

statusint
holds the entry status. Possible values are

1

an initial entry with u set to $b$.

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$.

powerfloat
holds the regularization power, $p \geq 2$.

weightfloat
holds the strictly positive regularization weight, $\sigma$.

undarray(m)
holds the result vector when initial or return status = 1, 2 or 4 (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 lsrt.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] lsrt.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, $p \geq 2$ or $\sigma > 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.

bitersint
the total number of inner iterations performed.

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

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

objfloat
the value of the objective function.

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

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.

lsrt.terminate()#

Deallocate all internal private storage.

example code#

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

# set problem data

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

# set default options
options = lsrt.initialize()

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

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

# problem

print(" solve problem")
status = 1
weight = 1.0
u = b
while status > 0:
  status, x, u, v = lsrt.solve_problem(status, m, n, power, weight, 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 = lsrt.information()
print(" optimal f: %.4f" % inform['obj'])
#print(" x:",x)
print('** lsrt exit status:', inform['status'])

# deallocate internal data

lsrt.terminate()

This example code is available in $GALAHAD/src/lsrt/Python/test_lsrt.py .