GALAHAD SHA package#

purpose#

The sha package finds a component-wise secant approximation to the Hessian matrix $H(x)$, for which $(H(x))_{i,j} = \partial f^2 (x) / \partial x_i \partial x_j$, $1 \leq i, j \leq n$, using values of the gradient $g(x) = \nabla_x f(x)$ of the function $f(x)$ of $n$ unknowns $x = (x_1, \ldots, x_n)^T$ at a sequence of given distinct $\{x^{(k)}\}$, $k \geq 0$. More specifically, given differences

\[s^{(k)} = x^{(k+1)} - x^{(k)} \;\;\mbox{and}\;\; y^{(k)} = g(x^{(k+1)}) - g(x^{(k)})\]

the package aims to find an approximation $B$ to $H(x)$ for which the secant conditions $B s^{(k)} \approx y^{(k)}$ hold for a chosen set of values $k$. The methods provided take advantage of the entries in the Hessian that are known to be zero.

The package is particularly intended to allow gradient-based optimization methods, that generate iterates $x^{(k+1)} = x^{(k)} + s^{(k)}$ based upon the values $g( x^{(k)})$ for $k \geq 0$, to build a suitable approximation to the Hessian $H(x^{(k+1)})$. This then gives the method an opportunity to accelerate the iteration using the Hessian approximation.

See Section 4 of $GALAHAD/doc/sha.pdf for additional details.

method#

The package computes the entries in the each row of $B$ one at a time. The entries $b_{ij}$ in row $i$ may be chosen to

\[(1) \;\;\; \min_{b_{i,j}}\!\mbox{imize} \;\; \sum_{k \in {\cal I}_i} \left[ \sum_{{\scriptscriptstyle \mbox{nonzeros}}\; j} b_{i,j} s_j^{(k)} - y_i^{(k)} \right]^2,\]

where ${\cal I}_i$ is ideally chosen to be sufficiently large so that (1) has a unique minimizer. Since this requires that there are at least as many $(s^{(k)}, y^{(k)})$ pairs as the maximum number of nonzeros in any row, this may be prohibitive in some cases. We might then be content with a minimum-norm (under-determined) least-squares solution; each row may then be processed in parallel. Or, we may take advantage of the symmetry of the Hessian, and note that if we have already found the values in row $j$, then the value $b_{i,j} = b_{j,i}$ in (1) is known before we process row $i$. Thus by ordering the rows and exploiting symmetry we may reduce the numbers of unknowns in future unprocessed rows.

In the analysis phase, we order the rows by constructing the connectivity graph—a graph comprising nodes $1$ to $n$ and edges connecting nodes $i$ and $j$ if $h_{i,j}$ is everywhere nonzero—of $H(x)$. The nodes are ordered by increasing degree (that is, the number of edges emanating from the node) using a bucket sort. The row chosen to be ordered next corresponds to a node of minimum degree, the node is removed from the graph, the degrees updated efficiently, and the process repeated until all rows have been ordered. This often leads to a significant reduction in the numbers of unknown values in each row as it is processed in turn, but numerical rounding can lead to inaccurate values in some cases. A useful remedy is to process all rows for which there are sufficient $(s^{(k)}, y^{(k)})$ as before, and then process the remaining rows taking into account the symmetry. That is, the rows and columns are rearranged so that the matrix is in block form

\[\begin{split}B = \begin{pmatrix} B_{11} & B_{12} \\ B^T_{12} & B_{22}\end{pmatrix},\end{split}\]

the $( B_{11} \;\; B_{12})$ rows are processed without regard for symmetry but give the $2,1$ block $B^T_{12}$, and finally the $2,2$ block $B_{22}$ is processed knowing $B^T_{12}$ again without respecting symmetry. The rows in blocks $( B_{11} \;\; B_{12})$ and $B_{22}$ may be computed in parallel. It is also possible to generalise this so that $B$ is decomposed into $r$ blocks, and the blocks processed one at a time recursively using the symmetry from previos rows. More details of the precise algorithms (Algorithms 2.1–2.5) are given in the reference below. The linear least-squares problems (1) themselves are solved by a choice of LAPACK packages.

references#

The method is described in detail in

J. M. Fowkes, N. I. M. Gould and J. A. Scott, “Approximating large-scale Hessians using secant equations”. Preprint TR-2024-001, Rutherford Appleton Laboratory.

parametric real type T and integer type INT#

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). The symbol INT refers to a parametric integer type that may be Int32 (32-bit integer) or Int64 (64-bit integer).

callable functions#

    function sha_initialize(T, INT, 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 sha_control_type)

status

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

0
The initialization was successful.

    function sha_read_specfile(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/sha/SHA.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/sha.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 sha_control_type)
specfile	is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file

    function sha_analyse_matrix(T, INT, control, data, status, n, ne, row, col, m)

Analsyse the sparsity structure of $H$ to generate information that will be used when estimating its values.

Parameters:

control	is a structure whose members provide control paramters for the remaining prcedures (see sha_control_type)
data	holds private internal data
status	is a scalar variable of type INT that gives the exit status from the package. Possible values are: 0 The import and analysis were conducted successfully. 1 Insufficient data pairs $(s_i,y_i)$ have been provided, as m is too small. The returned $B$ is likely not fully accurate. -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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 A restriction n > 0, ne $\geq$ 0 or 0 $\leq$ row[i] $\leq$ col[i] $\leq$ n has been violated.
n	is a scalar variable of type INT, that holds the number of variables
ne	is a scalar variable of type INT, that holds the number of entries in the upper triangular part of $H$.
row	is a one-dimensional array of size ne and type INT, that holds the row indices of the upper triangular part of $H$.
col	is a one-dimensional array of size ne and type INT, that holds the column indices of the upper triangular part of $H$.
m	is a scalar variable of type INT, that gives the minimum number of $(s^{(k)},y^{(k)})$ pairs that will be needed to recover a good Hessian approximation.

    function sha_recover_matrix(T, INT, data, status, ne, m, ls1, ls2, strans,
                                ly1, ly2, ytrans, val, order)

Estimate the nonzero entries of the Hessian $H$ by component-wise secant approximation.

Parameters:

data	holds private internal data
status	is a scalar variable of type INT that gives the exit status from the package. Possible values are: 0 The recovery was successful. -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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. -31 sha.recover_matrix has been called before sha.analyse_matrix.
ne	is a scalar variable of type INT that holds the number of entries in the lower triangular part of the symmetric matrix $A$.
val	is a one-dimensional array of size ne and type T that holds the values of the entries of the lower triangular part of the symmetric matrix $A$ in any of the supported storage schemes.
ne	is a scalar variable of type INT, that holds the number of entries in the upper triangular part of $H$.
m_available	is a scalar variable of type INT, that holds the number of differences provided. Ideally this will be as large as m as reported by sha_analyse_matrix, but better still there should be a further control.extra_differences to allow for unlikely singularities.
ls1	is a scalar variable of type INT, that holds the leading (first) dimension of the array strans.
ls2	is a scalar variable of type INT, that holds the trailing (second) dimension of the array strans.
strans	is a two-dimensional array of size [ls1][ls2] and type T, that holds the values of the vectors $\{s^{(k) T}\}$. Component [$k$][$i$] should hold $s_i^{(k)}$.
ly1	is a scalar variable of type INT, that holds the leading (first) dimension of the array ytrans.
ly2	is a scalar variable of type INT, that holds the trailing (second) dimension of the array ytrans.
ytrans	is a two-dimensional array of size [ly1][ly2] and type T, that holds the values of the vectors $\{y^{(k) T}\}$. Component [$k$][$i$] should hold $y_i^{(k)}$.
val	is a one-dimensional array of size ne and type T, that holds the values of the entries of the upper triangular part of the symmetric matrix $H$ in the sparse coordinate scheme.
order	is a one-dimensional array of size m and type INT, that holds the preferred order of access for the pairs $\{(s^{(k)},y^{(k)})\}$. The $k$-th component of order specifies the row number of strans and ytrans that will be used as the $k$-th most favoured. order need not be set if the natural order, $k, k = 1,...,$ m, is desired, and this case order should be C_NULL.

    function sha_information(T, INT, data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see sha_inform_type)

status

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

0
The values were recorded successfully

    function sha_terminate(T, INT, data, control, inform)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a structure containing control information (see sha_control_type)
inform	is a structure containing output information (see sha_inform_type)

available structures#

sha_control_type structure#

    struct sha_control_type{INT}
      f_indexing::Bool
      error::INT
      out::INT
      print_level::INT
      approximation_algorithm::INT
      dense_linear_solver::INT
      extra_differences::INT
      sparse_row::INT
      recursion_max::INT
      recursion_entries_required::INT
      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

INT error

error and warning diagnostics occur on stream error

INT out

general output occurs on stream out

INT print_level

the level of output required. <= 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

INT approximation_algorithm

which approximation algorithm should be used?

1 : unsymmetric, parallel (Alg 2.1 in paper)
2 : symmetric (Alg 2.2 in pape)
3 : composite, parallel (Alg 2.3 in paper)
4 : composite, block parallel (Alg 2.4 in paper)

INT dense_linear_solver

which dense linear equation solver should be used?

1 : Gaussian elimination
2 : QR factorization
3 : singular-value decomposition
4 : singular-value decomposition with divide-and-conquer

INT extra_differences

if available use an addition extra_differences differences

INT sparse_row

a row is considered sparse if it has no more than .sparse_row entries

    ipc_ recursion_max

limit on the maximum number of levels of recursion (Alg. 2.4)

    ipc_ recursion_entries_required

the minimum number of entries in a reduced row that are required if a further level of recuresion is allowed (Alg. 2.4)

Bool space_critical

if space is critical, ensure allocated arrays are no bigger than needed

Bool deallocate_error_fatal

exit if any deallocation fails

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’

sha_inform_type structure#

    struct sha_inform_type{T,INT}
      status::INT
      alloc_status::INT
      max_degree::INT
      differences_needed::INT
      max_reduced_degree::INT
      approximation_algorith_used::INT
      bad_row::INT
      max_off_diagonal_difference::T
      bad_alloc::NTuple{81,Cchar}

detailed documentation#

inform derived type as a Julia structure

components#

INT status

return status. See SHA_solve for details

INT alloc_status

the status of the last attempted allocation/deallocation.

INT max_degree

the maximum degree in the adgacency graph.

INT differences_needed

the number of differences that will be needed.

INT max_reduced_degree

the maximum reduced degree in the adgacency graph.

INT approximation_algorithm_used

the approximation algorithm actually used

INT bad_row

a failure occured when forming the bad_row-th row (0 = no failure).

T max_off_diagonal_difference

the maximum difference between estimated Hessian off-diagonal pairs if approximation algorithm 1, 3 or 4 has been employed and control.average_off_diagonals is true. It will be zero otherwise.

NTuple{81,Cchar} bad_alloc

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

Index