NOTE 224 – AIPS++ Least Squares background

Wim Brouw

22 January 1999

A pdf version of this note is available.

1 Introduction
2 Linear equations
2.1 Real
  2.1.1 Errors
2.2 Complex
  2.2.1 Separable complex
  2.2.2 Errors
3 Dependant linear equations
3.1 Unknown constraints
  3.1.1 Errors
3.2 Known constraints
  3.2.1 Errors
4 Non-linear equations
4.1 Errors
5 Summary

1 Introduction

Trying to ﬁnd a more stable non-linear least squares method (LSQ) than the one currently used in Newstar, I found another few areas that could be added or improved in the existing routines. In the end I produced a set of routines to handle linear and non-linear cases, both with or without SVD, and constraint equations.

In addition the present document describes the usage and background of the diﬀerent routines.

Least squares solutions in Aips++ have the following general background:

all are based on creating normal equations from an, initially unknown, number of condition equations
the actual LSQ-object is a symmetric (or Hermitian) matrix
matrix inversions are done using in-place Cholesky methods where possible: it is for standard use the fastest, least memory consuming and a stable method. In cases where the inverse matrix is needed, the method uses resources of the same order as full LU-decomposition (only if explicit constraint equations are present)
input/output to the LSQ-object is done in any precision, internally all calculations are done in double precision
complex numbers are assumed to be contiguous pairs of real numbers with the real part being the ﬁrst

The following terminology is used throughout:

$n$: number of unknowns to be solved
$m$: number of simultaneous knowns
$N$: number of condition equations
$χ^{2}$: merit function
$z^{*}$: complex conjugate of $z$
$z_{ℜ}$: real part of $z$
$z_{ℑ}$: imaginary part of $z$
$x$: column vector of unknowns $x_{0}, \dots, x_{n - 1}$
$a_{i}$: vector of factors $a_{0, i}, \dots, a_{n - 1, i}$ in condition equation $i$
( $i = 0, \dots, N - 1$ )
$C$: the $N \times n$ array of condition equations
$A$: $n \times n$ array: normal equations matrix
$L$: $n$ column vector: right-hand side normal equations
$y_{i}$: model value for condition equation $i$ ( $i = 0, \dots, N - 1$ )
$l_{i}$: measured value for condition equation $i$ ( $i = 0, \dots, N - 1$ )
$σ_{i}$: standard deviation for condition equation $i$ ( $i = 0, \dots, N - 1$ )
$w_{i}$: weight for condition equation $i$ ( $w_{i} = 1 ∕ σ_{i}^{2}$ )
$[a b]$: shorthand for $\sum_{i = 0}^{N - 1} w_{i} a_{i} b_{i}$

The condition equations can, with the above deﬁnitions, be written in one of the following ways:

a_{i} \cdot x = l_{i} i = 0, \dots, N - 1

(1)

or as:

\begin{array}{rcl} \sum_{k = 0}^{n - 1} a_{i k} x_{k} = l_{i} i = 0, \dots, N - 1 & (2) \end{array}

or as:

C \cdot x = l

(3)

The normal matrix is deﬁned as:

A = C^{T} \cdot Q^{- 1} \cdot C

(4)

resulting in the normal equations:

A \cdot x = L

(5)

where $L = C^{T} \cdot Q^{- 1} \cdot l$ and $Q$ is the covariance matrix of the observations. Here only covariance matrices with the same value in each column $i$ (the ‘weight’ $w_{i}$ ) are considered.

2 Linear equations

2.1 Real

The merit function we want to minimise is:

\begin{array}{rcl} χ^{2} = \sum_{i = 0}^{N - 1} {[\frac{l_{i} - y_{i}}{σ_{i}}]}^{2} & (6) \end{array}

For the minimum of $χ^{2}$ holds:

\frac{\partial χ^{2}}{\partial x_{k}} = 0 k = 0, \dots, n - 1

(7)

which leads to the following set of normal equations:

\sum_{i = 0}^{N - 1} \frac{[l_{i} - y_{i}]}{σ_{i}^{2}} \frac{\partial y_{i}}{\partial x_{k}} = 0 k = 0, \dots, n - 1

(8)

or in matrix form:

(\begin{matrix} [a_{0} a_{0}] & [a_{0} a_{1}] & \dots & [a_{0} a_{n - 1}] \\ [a_{1} a_{0}] & [a_{1} a_{1}] & \dots & [a_{1} a_{n - 1}] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ [a_{n - 1} a_{0}] & [a_{n - 1} a_{1}] & \dots & [a_{n - 1} a_{n - 1}] \end{matrix}) \cdot x = (\begin{matrix} [a_{0} l] \\ [a_{1} l] \\ ⋮ \\ [a_{n - 1} l] \end{matrix})

(9)

or:

A \cdot x = L

(10)

This symmetric set of equations can be solved, if the matrix is positive deﬁnite, i.e. if there are no dependencies between the columns. The solution is given by:

x = A^{- 1} \cdot L

(11)

2.1.1 Errors

After solution for the unknowns $x$ , $χ^{2}$ as deﬁned in (6) can be directly calculated if we rewrite it as:

χ^{2} = [l l] - 2 \sum_{k = 0}^{n - 1} [a_{k} l] x_{k} + \sum_{i = 0}^{N - 1} w_{i} y_{i}^{2}

(12)

Noting that $y_{i} = \sum_{k = 0}^{n - 1} a_{k} x_{k}$ and using equation (9) to note that $[a_{i} l] = \sum_{k = 0}^{n - 1} [a_{i} a_{k}] x_{k}$ , we can rewrite (12) as:

χ^{2} = [l l] - \sum_{k = 0}^{n - 1} x_{k} [a_{k} l]

(13)

The $χ^{2}$ could be used to assess the goodness of ﬁt if the actual $σ_{i}$ ’s were known, and the errors are normal distributed. In general the actual values of $σ_{i}$ are not known, and often the distribution is not normal (.e.g. if we solve for logarithmic values). An estimate of the standard deviation can be made by:

σ^{2} = \frac{[{(l_{i} - y_{i})}^{2}]}{N - n}

(14)

which can be estimated by:

σ_{o}^{2} = \frac{χ^{2}}{N - n}

(15)

to give ‘an error per observation’. The ‘error per unit weight’, or the standard deviation, can be expressed as:

σ_{w}^{2} = \frac{χ^{2}}{[1]} \frac{N}{N - n}

(16)

The uncertainty in the solution $x_{i}$ can be expressed as:

σ^{2} (x_{i}) = \sum_{k = 1}^{N - 1} σ_{k}^{2} {(\frac{\partial x_{i}}{\partial y_{k}})}^{2}

(17)

Since

x_{i} = \sum_{k = 0}^{n - 1} {(A^{- 1})}_{i k} [a_{k} l]

(18)

we have:

\frac{\partial x_{i}}{\partial y_{k}} = \sum_{j = 0}^{n - 1} {(A^{- 1})}_{i j} \frac{a_{k j}}{σ_{k}^{2}}

(19)

leading to:

σ^{2} (x_{i}) = \sum_{k = 0}^{n - 1} \sum_{l = 0}^{n - 1} {(A^{- 1})}_{i k} {(A^{- 1})}_{i l} [\sum_{j = 0}^{N - 1} \frac{a_{k j} a_{l j}}{σ_{j}^{2}}]

(20)

Doing the sums, this equation reduces to:

σ^{2} (x_{i}) = {(A^{- 1})}_{i i}

(21)

If the $σ_{i}$ was not known originally, the estimate of the standard uncertainties in the unknowns $x$ are:

σ (x_{i}) = σ_{o} \sqrt{{(A^{- 1})}_{i i}}

(22)

If the uncertainties in the unknowns are wanted, the inverse matrix $A^{- 1}$ is calculated by backsubstitution of the identity matrix, and the uncertainties by (22).

2.2 Complex

The merit function we want to minimise in this case is:

χ^{2} = \sum_{i = 0}^{N - 1} [\frac{l_{i} - y_{i}}{σ_{i}}] {[\frac{l_{i} - y_{i}}{σ_{i}}]}^{*}

(23)

Diﬀerentiating $χ^{2}$ leads to the following set of normal equations:

(\begin{matrix} [a_{0}^{*} a_{0}] & [a_{0}^{*} a_{1}] & \dots & [a_{0}^{*} a_{n - 1}] \\ [a_{1}^{*} a_{0}] & [a_{1}^{*} a_{1}] & \dots & [a_{1}^{*} a_{n - 1}] \\ ⋮ & ⋮ & ⋱ & ⋮ \\ [a_{n - 1}^{*} a_{0}] & [a_{n - 1}^{*} a_{1}] & \dots & [a_{n - 1}^{*} a_{n - 1}] \end{matrix}) \cdot x = (\begin{matrix} [a_{0}^{*} l] \\ [a_{1}^{*} l] \\ ⋮ \\ [a_{n - 1}^{*} l] \end{matrix})

(24)

The normal matrix is Hermitian. It can be solved by Choleski methods. However, internally the matrix is converted to a real form. Although this has an, in general, small memory penalty, it has no inﬂuence on CPU time, and makes it possible to use the same routines for complex and real solutions. The conversion to real is done by splitting each element $A_{i j}$ of the normal matrix into ${A_{ℜ}}_{i j} + ı {A_{ℑ}}_{i j}$ and replacing it by:

A_{i j} = (\begin{matrix} {A_{i j}}_{ℜ} & - {A_{i j}}_{ℑ} \\ {A_{i j}}_{ℑ} & {A_{i j}}_{ℜ} \end{matrix})

(25)

and simular replacements for the vector elements $x_{i}$ and $L_{i}$ as:

x_{i} = (\begin{matrix} {x_{i}}_{ℜ} \\ {x_{i}}_{ℑ} \end{matrix})

(26)

and:

L_{i} = (\begin{matrix} {L_{i}}_{ℜ} \\ {L_{i}}_{ℑ} \end{matrix})

(27)

Another reason for solving real rather than complex equations is given in the next section.

2.2.1 Separable complex

In cases where both the unknowns $x$ and their complex conjugates $x^{*}$ appear in the condition equations, diﬀerentiating the merit function (23) will not produce a symmetric or Hermitian normal matrix, since there exists no linear relation between $x_{i}$ and $x_{i}^{*}$ . We could, of course, solve for $2 n$ complex equations with added constraints that the sum of even and odd unknowns must be real, and their diﬀerence imaginary, but this will lead to $4 n$ complex equations.

If, however, we consider each complex unknown as two real unknowns (i.e. ${x_{i}}_{ℜ}$ and ${x_{i}}_{ℑ}$ ) then diﬀerentiating (23) produces the following symmetric set of $2 n$ real equations:

(\begin{matrix} {[{a_{0}}^{*} a_{0}]}_{ℜ} & {[{a_{1}}^{*} a_{0}]}_{ℑ} & {[{a_{2}}^{*} a_{0}]}_{ℜ} & \dots & {[{a_{2 n - 1}}^{*} a_{0}]}_{ℑ} \\ - {[{a_{0}}^{*} a_{1}]}_{ℑ} & {[{a_{1}}^{*} a_{1}]}_{ℜ} & - {[{a_{2}}^{*} a_{1}]}_{ℑ} & \dots & {[{a_{2 n - 1}}^{*} a_{1}]}_{ℜ} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - {[{a_{0}}^{*} a_{2 n - 1}]}_{ℑ} & {[{a_{1}}^{*} a_{2 n - 1}]}_{ℜ} & - {[{a_{2}}^{*} a_{2 n - 1}]}_{ℑ} & \dots & {[{a_{2 n - 1}}^{*} a_{2 n - 1}]}_{ℜ} \end{matrix}) \cdot

\cdot (\begin{matrix} {x_{0}}_{ℜ} \\ {x_{0}}_{ℑ} \\ {x_{1}}_{ℜ} \\ ⋮ \\ {x_{n - 1}}_{ℑ} \end{matrix}) = (\begin{matrix} {[{a_{0}}^{*} l]}_{ℜ} \\ {[{a_{1}}^{*} l]}_{ℑ} \\ {[{a_{2}}^{*} l]}_{ℜ} \\ ⋮ \\ {[{a_{2 n - 1}}^{*} l]}_{ℑ} \end{matrix})

(28)

A number of special cases can be distinguished.

In the ‘normal’ complex case (previous section), $a_{2 i} \equiv a_{2 i + 1}$ , and (28) reduces to (25).

If all $a_{2 i}$ and $a_{2 i + 1}$ are real, but speciﬁed, e.g., as a complex number, (28) deteriorates into:

(\begin{matrix} [{a_{0}}_{ℜ} {a_{0}}_{ℜ}] & 0 & [{a_{0}}_{ℜ} {a_{1}}_{ℜ}] & \dots & 0 \\ 0 & [{a_{0}}_{ℑ} {a_{0}}_{ℑ}] & 0 & \dots & [{a_{0}}_{ℑ} {a_{n - 1}}_{ℑ}] \\ [{a_{1}}_{ℜ} {a_{0}}_{ℜ}] & 0 & [{a_{1}}_{ℜ} {a_{1}}_{ℜ}] & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & [{a_{n - 1}}_{ℑ} {a_{0}}_{ℑ}] & 0 & \dots & [{a_{n - 1}}_{ℑ} {a_{n - 1}}_{ℑ}] \end{matrix}) \cdot

\cdot (\begin{matrix} {x_{0}}_{ℜ} \\ {x_{0}}_{ℑ} \\ {x_{1}}_{ℜ} \\ ⋮ \\ {x_{n - 1}}_{ℑ} \end{matrix}) = (\begin{matrix} [{a_{0}}_{ℜ} l_{ℜ}] \\ [{a_{0}}_{ℑ} l_{ℑ}] \\ [{a_{1}}_{ℜ} l_{ℜ}] \\ ⋮ \\ [{a_{n - 1}}_{ℑ} l_{ℑ}] \end{matrix})

(29)

Note that in this case there is a true separation of the real and imaginary parts of the diﬀerent unknowns, and two separate real solutions with each $n$ unknowns will produce exactly the same result.

The full and special cases are all catered for in the aips++ routines.

2.2.2 Errors

Using arguments comparable to those in (2.1.1), we can get (13) for the complex case as:

χ^{2} = [l l^{*}] - \sum_{k = 0}^{n - 1} ({x_{k}}_{ℜ} {[a_{2 k}^{*} l]}_{ℜ} + {x_{k}}_{ℑ} {[a_{2 k + 1}^{*} l]}_{ℑ})

(30)

with the $a$ ’s and $x$ ’s as deﬁned in (29).

3 Dependant linear equations

If there are not enough independent condition equations, the normal matrix $A$ cannot be inverted, and a call to WNMLTN will fail with a .false. return value.

The equations could still be solved if some additional ‘constraint’ equations would be introduced. In the more complex cases the precise, let alone the best, form for these additional equations is diﬃcult to determine (e.g. the redundancy situation in Westerbork).

A method known as ‘Singular value decomposition’ (SVD) can be used to obtain the minimal set of orthogonal equations that have to be added to solve the LSQ problem. Several implementations exist in the literature.

In general we can distinguish three types of constrained equations:

the minimum number of suﬃcient constraint equations are known to be able to solve a system of equations
constraints are used to add additional information (e.g. the sum of angles of a triangle)
no actual constraint equations are known

All three cases are handled in the LSQ package, the ﬁrst two in the same way.

The general constraint situation arises from the use of Lagrange multiplicators. Assume that in addition to the condition equations, with measured values, we have a set of $p$ rigorous equations:

ϕ_{i} (x) = 0 i = 0, \dots, p - 1

(31)

We must therefore make $χ^{2}$ minimal, subject to the set of (31), or:

\sum_{i = 0}^{n - 1} \frac{\partial χ^{2}}{\partial x_{i}} d x_{i} = 0

(32)

subject to the conditions:

\sum_{i = 0}^{n - 1} \frac{\partial ϕ_{k}}{\partial x_{i}} d x_{i} = 0 k = 0, \dots, p - 1

(33)

which leads to a set of $n + p$ equations:

\frac{\partial χ^{2}}{\partial x_{i}} + \sum_{k = 0}^{p - 1} λ_{k} \frac{\partial ϕ_{k}}{\partial x_{i}} = 0 i = 0, \dots, n - 1

(34)

together with the (31).

Note that I have chosen for having constraint equations linear in the unknowns, with a zero value. In cases where this is not adequate (e.g. $x + y + z = 360$ ) a simple linear transformation will suﬃce to make it e.g. $x^{'} + y^{'} + z^{'} = 0$ .

Deﬁning the second term in (34) as $B_{i k}$ , we can write our expanded set of normal equations as:

(\begin{matrix} A & B \\ B^{T} & 0 \end{matrix}) (\begin{matrix} x \\ λ \end{matrix}) = (\begin{matrix} L \\ 0 \end{matrix})

(35)

3.1 Unknown constraints

The routines to handle this situation use a method based on “Pseudo-inverse berekening en Cholesky factorisatie”, Hans van der Marel, 27 maart 1990, Faculty of Geodesy, Delft University. I will brieﬂy describe Hans’ paper, together with the additions I have made.

In the case we are considering, the normal equation $A$ has not full column rank. Therefore, there exists, if the rank defect is $p$ , an $n \times p$ matrix $G$ , a basis for the null-space of $A$ , with $A G = 0$ . If we assume that $B$ has just suﬃcient constraint equations to solve the rank defect, the inverse of the matrix in (35) is:

{(\begin{matrix} A & B \\ B^{T} & 0 \end{matrix})}^{- 1} = (\begin{matrix} A^{#} & G {(B^{T} G)}^{- 1} \\ {(G^{T} B)}^{- 1} G^{T} & 0 \end{matrix})

(36)

with $A^{#}$ the pseudo-inverse of $A$ . A number of relations hold for $A^{#}$ , the most important for us that $B^{T} A^{#} = 0$ and $A^{#} B = 0$ , or $B$ is a base for the null-space of $A^{#}$ . For a mininorm solution we can choose $B = G$ . In that case $A^{#}$ is the Moore-Penrose inverse, and $B^{T} G$ is regular and symmetric.

We can now proceed in the following way:

Do a Cholesky factorisation of the normal equations

A

until the remaining columns of

A

are dependent. Pivoting along the diagonal of

A

occurs to make sure that

A

can be partitioned in an independent and a dependent part with rank defect

n_{2}

. Dependency is determined by looking at the angle between a column and the space formed by the already determined independent space (the so-called ‘collinearity number’).

If $n_{1} = n - n_{2}$ , we can say that after $n_{1}$ Choleski steps, we have:

A = (\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}) = (\begin{matrix} U_{11}^{T} & 0 \\ U_{12}^{T} & I \end{matrix}) \cdot (\begin{matrix} U_{11} & U_{12} \\ 0 & {\bar{A}}_{22} \end{matrix})

(37)

with $U_{11}$ the Cholesky factor of $A_{11}$ , $U_{12} = {U_{11}^{T}}^{- 1} A_{12}$ and ${\bar{A}}_{22} = A_{22} - U_{12}^{T} U_{12}$ .

The collinearity angle $δ$ can be determined for the current column $i$ by: ${sin}^{2} δ = u_{i i}^{2} ∕ a_{i i}$ .

Determine a $G_{1}$ replacing the (rectangular) $U_{12}$ from $U_{11} G_{1} = - U_{12}$ .
Determine a (symmetric) $G_{2}$ replacing $A_{22}$ from the rank basis $I + G_{1}^{T} G_{1} = (G_{2}^{T} - 1) (G_{2} - 1)$
determine constraint equations $G = {(G_{1} G_{2}, G_{2})}^{T}$ .
Solve for the $n_{1}$ independent unknowns $x_{1}$ using $A_{11}$ and $L_{1}$ by back substitution
Solve constraint equations for the remaining $n_{2}$ unknowns, using $G_{2}$ by back substitution: $x_{2} = - G_{2}^{- 1} G_{1}^{T} x_{1}$
Make Baarda’s S-transform of independent $n_{1}$ unknowns: $x_{1} = x_{1} + G_{1} x_{2}$ .

Setting $F = A_{1}^{- 1} L_{1}$ , $D = - G_{2}^{- 1} G_{1}^{T}$ and $E = (\begin{matrix} I_{n_{1}} + G_{1} D \\ D \end{matrix})$ , we can write the above steps as:

(\begin{matrix} x_{1} \\ x_{2} \end{matrix}) = (\begin{matrix} I_{n_{1}} + G_{1} D \\ D \end{matrix}) (\begin{matrix} A_{1}^{- 1} L_{1} \end{matrix}) = E \cdot F

(38)

3.1.1 Errors

Errors are determined in the same way as described in (2.1.1), where it should be noted that (13) can either be used summing over $n_{1}$ variables and using the ﬁrst guess of $x_{1}$ , or over $n$ and using the ﬁnal $x$ . Internally the ﬁrst option is used. The uncertainties in $x$ are determined by calculating the covariance matrix. Using similar arithmetic as in (2.1.1), it can be shown that:

A^{- 1} = E \cdot a_{1}^{- 1} \cdot E^{T}

(39)

$A^{- 1}$ is calculated by ﬁrst solving $H = E (A_{1}^{- 1} \cdot I_{n_{1}})$ and then $A^{- 1} = E \cdot H^{T}$ .

3.2 Known constraints

In the case of known constraints, hence when $B$ is known in (35), this equation can be solved. Cholesky decomposition does not work in this case, and a, transparent, Crout LU decomposition is done to determine the (symmetric) covariance (i.e. inverse) matrix of the left hand side of (35). .

3.2.1 Errors

Errors are determined as explained in (2.1.1). Since all constraint equations are $\equiv 0$ , the sum in (13) is taken over $n$ , not over $n + p$ if $p$ are the number of constraint equations.

4 Non-linear equations

In the case of non-linear condition equations, no simple solution to minimise the merit function (e.g. (6)) exists. A simple solution is to take the condition equation, make a ﬁrst order Taylor expansion around an estimated $\hat{x}$ ; solve for $d x$ , and iterate till solution found.

However, better and more stable methods exist (e.g. steepest-descend) for some circumstances. A method that seems to be quite stable and using both steepest-descent and Taylor expansion, is the method by Levenberg-Marquardt (see e.g. “Numerical recipes”, W.H. Press et al., Cambridge University Press).

If we have an estimate for $x$ , we can ﬁnd a better one by:

{\hat{x}}_{n e x t} = \hat{x} + H^{- 1} \cdot [- \nabla χ^{2} (\hat{x})]

(40)

where $H$ is the Hessian matrix of $χ^{2}$ . If our approximation is not good enough, we could use the steepest-descent by calculating:

{\hat{x}}_{n e x t} = \hat{x} - c o n s t a n t \cdot [- \nabla χ^{2} (\hat{x})]

(41)

The Hessian matrix $H$ has as elements:

H_{i j} = \frac{\partial^{2} χ^{2}}{\partial x_{i} \partial x_{j}} = 2 \sum_{k = 0}^{N - 1} \frac{1}{σ_{i}^{2}} [\frac{\partial y_{k}}{\partial x_{i}} \frac{\partial y_{k}}{\partial x_{j}} - (l_{i} - y_{i}) \frac{\partial^{2} y_{i}}{\partial x_{i} \partial x_{j}}]

(42)

By dropping the second term in (42), and multiplying each $H_{i i}$ term with $(1 + λ)$ , and deﬁning the ﬁrst derivative of $χ^{2}$ as $b$ , we can combine the equations (40), (41) into:

H \cdot d x = b

(43)

which can be solved as standard normal equations.

Choosing $λ$ is the crux of the matter, and where to stop iterating. A start value of $λ = 0.001$ is used in the following routines. The method looks at the $χ^{2}$ for $\hat{x} + d x$ . If it has increased over the value of $χ^{2}$ for $\hat{x}$ , $\hat{x}$ is unchanged, and $λ = 10 λ$ . If there is a decrease; a new value for $\hat{x}$ is used, and $λ = λ ∕ 10$ . Iteration can be stopped if the fractional decrease in $χ^{2}$ is small (say $< 0.001$ ); never if there is an increase in $χ^{2}$ .

4.1 Errors

Errors are calculated as described in (2.1.1). The errors returned by WNMLSN are, of course, only approximations, since the original equations are non-linear, but give a good impression of the residuals. The covariance matrix should be calculated by doing a ﬁnal linear run on the residuals, and solve the equations.

5 Summary

An LSQ-object takes a total space of:

\begin{array}{rcl} 9 & + & (a d m i n i s t r a t i o n) \\ (n + p) (n + p + 1) ∕ 2 & + & (n o r m a l a r r a y) \\ 4 m & + & (e r r o r c a l c u l a t i o n s) \\ m (m + p) & + & (k n o w n s i d e s) \\ (n + p + 1) ∕ 2 & + & (p i v o t a r e a) \\ (n + p) & (s o l u t i o n) 8 b y t e w o r d s & (44) \end{array}

where $n$ is the number of unknowns ( $2 n$ if complex); $m$ the number of simultaneously to be solved equations; $p$ the number of external constraint equations. In the non-linear case two LSQ-objects are used. In cases where the inverted normal array is calculated (known constraints) a temporary storage of $n^{2}$ 8-byte words is used.

The following calls are available:

To be given

wnb, August 5, 2015