System of nonlinear equations with separable variables

Let f : R_{n} → R_{n} is a function of variable x ∈ R_{n}. Then, if α is a solution of the equation f(α) = 0, the function \frac{1}{2}∥f(x)∥^{2} = \frac{1}{2}f^{T}(x)f(x) has global minimum at the point α.

Let f : R_{n} → R_{m}, m > n, is a function of variable x ∈ R_{n}. Then, if the function \frac{1}{2}∥f(x)∥^{2} has global minimum at the point α, we say, that α is a solution of the equation f(α) = 0 according to the criterion of least sum of squares of residuals.

Let f : R_{n} → R_{m}, m ≥ n, is a function of variable x ∈ R_{n} such that

f(x) ≡ f( \begin{matrix}x_{1} \\ x_{2} \end{matrix} ) = b(x_{2}) − A(x_{2})x_{1} ,

x_{1} ∈ R_{k}, x_{2} ∈ R_{n−k}, A : R_{n−k} → R_{m×k}, b : R_{n−k} → R_{m} (the variables of the function f(x) are separable, i.e. f(x) is a linear function of x_{1} and nonlinear function of x_{2}). Let f_{2} : R_{n−k} → R_{m} is a function of x_{2} such that¶

f_{2}(x_{2}) = b(x_{2}) − A(x_{2})A^{+}(x_{2})b(x_{2}) = [1_{m} − A(x_{2})A^{+}(x_{2})] b(x_{2}) .

Then^{1}

f(α) ≡ f( \begin{matrix}α_{1} \\ α_{2} \end{matrix} ) = 0 ⇔ f_{2}(α_{2}) = 0 ∧ α_{1} = A^{+}(α_{2})b(α_{2}) .

In words, if α is a solution of the equation f(α) = 0†, α_{2} is a solution of the equation f_{2}(α_{2}) = 0 and α_{1} = A^{+}(α_{2})b(α_{2}). And also conversely (which is much more valuable from the practical point of view), if α_{2} is a solution of the equation f_{2}(α_{2}) = 0 and α_{1} = A^{+}(α_{2})b(α_{2}), α is a solution of the equation f(α) = 0§.

This theorem says, that if the variables of the function f(x) are separable, we can solve the equation f_{2}(α_{2}) = 0_{n−k} instead of the equation f(α) = 0_{n} and compute α_{1} “later”.

Let us suppose, that we use a derivative algorithm to solve the equation, which requires the knowledge of the Jacobi matrix J_{2}(x_{2}) ∼ (m, n − k) of the function f_{2}(x_{2})‡

J_{2}(x_{2}) = \{\frac{∂f_{2}_{i}}{∂x_{2}_{j}} \} , i = 1, …, m , j = 1, …, n − k

and also the gradient (i.e. the vector of the first derivatives) of the function \frac{1}{2}∥f_{2}(x_{2})∥^{2}

grad [\frac{1}{2}f_{2}^{T}(x_{2})f_{2}(x_{2})] = J_{2}^{T}(x_{2})f_{2}(x_{2}) .

However, usually we know nothing but the Jacobi matrix J(x) ∼ (m, n) of the function f(x)

J(x) = \{\frac{∂f_{i}}{∂x_{j}} \} , i = 1, …, m , j = 1, …, n ,

whose “left” submatrix of type (m, k) is the matrix A(x_{2}). The composition of the matrix J(x) and vector x_{1} = A^{+}(x_{2})b(x_{2}) leads to the matrix ( A(x_{2}) J_{n−k}(x_{2}) ).

¶ If f(x) = b − Ax_{1} = 0, then x_{1} = A^{+}b. If we introduce this solution “back” into the relation for f(x), we get f(x) = b − AA^{+}b = (1_{m} − AA^{+})b. In different words, the function f_{2}(x_{2}) is a composition of the function f(x) and vector x_{1} = A^{+}(x_{2})b(x_{2}).

† according to the criterion of least sum of squares of residuals

§ If the function \frac{1}{2}∥f(x)∥^{2} has global minimum at the point α, the function \frac{1}{2}∥f_{2}(x_{2})∥^{2} has global minimum at the point α_{2} and α_{1} = A^{+}(α_{2})b(α_{2}). If the function \frac{1}{2}∥f_{2}(x_{2})∥^{2} has global minimum at the point α_{2} and α_{1} = A^{+}(α_{2})b(α_{2}), the function \frac{1}{2}∥f(x)∥^{2} has global minimum at the point α.

‡ and that both the matrix A(x_{2}) and the vector b(x_{2}) are continuously differentiable at the point x_{2} and the matrix A(x_{2}) has continuous (i.e. actually constant, but \underline{not necessarily full}) rank r(A) ≤ n at the point x_{2}, thanks to which the matrix A^{+}(x_{2}) is continuously differentiable at the point x_{2}, too

In the following paragraph we will build a mathematical formalism, with the aid of which we will express the Jacobi matrix J_{2}(x_{2}) \underline{analytically} as a function of the matrix (A(x_{2}) J_{n−k}(x_{2}) ). However, first we will extend the validity of the multiplication and transposition operation of matrices to the “three-dimensional” case. Thus, let the multiplication of a three-dimensional and two-dimensional matrix takes its course in the first and second dimension, while the result is always a three-dimensional matrix (the second dimension can be equal to one). Analogically, let the transposition of a three-dimensional matrix consists in the exchange of the first and second dimension (the first and second index of all elements)¶.

Let f : R_{p} → R_{m×n} is a function of variable x ∈ R_{p} and D ≡ D_{x∈R_{p}} : R_{m×n} → R_{m×n×p} is a projection such that

D_{x∈R_{p}} [f(x)] = \{\frac{∂f_{ij}}{∂x_{k}} \} , i = 1, …, m , j = 1, …, n , k = 1, …, p .

Let A, B : R_{p} → R_{m×n}, a, b : R_{p} → R_{m×1} (i.e. a^{T}, b^{T} : R_{p} → R_{1×m}), a, b : R_{p} → R_{1×1} are functions of variable x ∈ R_{p}. Then

D(AB) = D(A)B + AD(B) ,
D(Ab) = D(A)b + AD(b) ,
D(a^{T}b) = D (a^{T}) a + a^{T}D(b) ,
D(aB) = D(a)B + aD(B) ,
D(ab) = D(a)b + aD(b) ,
D(ab) = D(a)b + aD(b) ,
D (A^{T}) = [D(A)]^{T} ,
D (a^{T}) = [D(a)]^{T} .

It is necessary to take into consideration that while the multiplication and transposition of matrices takes its course in the “paper plane”, the differentiation with respect to the elements of a vector goes in the third dimension, i.e. “perpendicularly to the paper”. Thus, the Jacobi matrix J(x) = D_{x∈R_{n}} [f(x)] = ( D_{x_{1}∈R_{k}} [f(x)] D_{x_{2}∈R_{n−k}} [f(x)] ) = ( A(x_{2}) D_{x_{2}∈R_{n−k}} [f(x)] ) is a matrix of type (m, 1, n), the Jacobi matrix J_{2}(x_{2}) = D_{x_{2}∈R_{n−k}} [f_{2}(x_{2})] is a matrix of type (m, 1, n − k) and the gradient of the function \frac{1}{2}f_{2}^{T}(x_{2})f_{2}(x_{2}) is a vector of type (1, 1, n − k) in this and several following paragraphs. Obviously

D_{x_{2}∈R_{n−k}} [f(x)] = D(b) − D(A)x_{1} .

The composition of this matrix† and vector x_{1} = A^{+}(x_{2})b(x_{2}) leads to the matrix J_{n−k}(x_{2}). Then§

J_{n−k}(x_{2}) = D(b) − D(A)A^{+}b ,

J_{2}(x_{2}) = D [f_{2}(x_{2})] = D [(1_{m} − AA^{+})b] = (1_{m} − AA^{+})D(b) − [D(A)A^{+} + AD(A^{+})] b =

= D(b) − D(A)A^{+}b − AA^{+}D(b) − AD(A^{+})b = J_{n−k}(x_{2}) − A [A^{+}D(b) + D(A^{+})b] .

¶ AB = C, Ab = c, A ∼ (m, n, p), B ∼ (n, k), C ∼ (m, k, p), b ∼ (n, 1), c ∼ (m, 1, p), C^{T} ∼ (k, m, p), c^{T} ∼ (1, m, p). Thus, in the first example (AB = C) we multiply sequentially p matrices of type (m, n) with one and the same matrix B of type (n, k) obtaining p matrices of type (m, k).

† I will say once more, that it is the “right” submatrix of type (m, 1, n − k) of the Jacobi matrix J(x)

§ hitherto we had to distinguish between D_{x∈R_{n}}, D_{x_{1}∈R_{k}} and D_{x_{2}∈R_{n−k}}, starting from this point we solely differentiate with respect to x_{2} ∈ R_{n−k}

From the definition of the Moore-Penrose pseudoinverse matrix¶ if follows that

A^{T}f_{2}(x_{2}) = A^{T}(1_{m} − AA^{+})b = [A^{T} − A^{T} (AA^{+})^{T}] b = [A^{T} − (AA^{+}A)^{T}] b = 0

and the gradient of the function \frac{1}{2}f_{2}^{T}(x_{2})f_{2}(x_{2}) can easily be obtained by mere “replacement” of the Jacobi matrix J_{2}(x_{2}) with the already well-known matrix J_{n−k}(x_{2})

grad [\frac{1}{2}f_{2}^{T}(x_{2})f_{2}(x_{2})] = J_{n−k}^{T}(x_{2})f_{2}(x_{2}) − [A^{+}D(b) + D(A^{+})b]^{T} A^{T}f_{2}(x_{2}) = J_{n−k}^{T}(x_{2})f_{2}(x_{2}) .

No relation between the matrices J_{2}(x_{2}) and J_{n−k}(x_{2}) follows from this fact yet. However, it is clear that both the shift vector δ_{2} = −J_{2}^{+}(x_{2})f_{2}(x_{2}) (if m = n, we talk about the Newton- -Raphson method, if m > n, we talk about the Gauss-Newton method) and the shift vector δ_{2} = −J_{n−k}^{+}(x_{2})f_{2}(x_{2}) are descent vectors of the function \frac{1}{2}∥f_{2}(x_{2})∥^{2} at the point x_{2} ≠ α_{2}, i.e.

\{grad [\frac{1}{2}f_{2}^{T}(x_{2})f_{2}(x_{2})]\}^{T} δ_{2} < 0 ,

if the matrices J_{2}(x_{2}) and J_{n−k}(x_{2}) have full rank† at the point x_{2}. Thus, if we substitute the matrix J_{n−k}(x_{2}) for the Jacobi matrix J_{2}(x_{2}) in the Newton-Raphson or Gauss-Newton method, both the methods will still converge to the right solution. But usually much more slowly (even in comparison with methods, which adopt the shift vector δ = −J^{+}(x)f(x) to solve the equation f(α) = 0, if we calculate the initial estimate of α_{1} from the initial estimate of α_{2} according to the relation x_{1} = A^{+}(x_{2})b(x_{2})), thanks to which the described simplification is almost unusable in practice§.

The relation between the matrices J_{2}(x_{2}) and J_{n−k}(x_{2}) can be derived even without an explicit expression of the three-dimensional matrix D(A^{+}) of type (k, m, n − k). It holds that

J_{2}(x_{2}) = D [f_{2}(x_{2})] = D [(1_{m} − AA^{+})b] = −D(AA^{+})b + (1_{m} − AA^{+})D(b)

and from the definition of the Moore-Penrose pseudoinverse matrix¶ it follows that^{1}

D(AA^{+}A) = D(A) = D(AA^{+})A + AA^{+}D(A) ,
D(AA^{+})AA^{+} = (1_{m} − AA^{+})D(A)A^{+} ,
D(AA^{+}AA^{+}) = D(AA^{+}) = D(AA^{+})AA^{+} + AA^{+}D(AA^{+}) ,

AA^{+}D(AA^{+}) = {[D(AA^{+})]^{T} (AA^{+})^{T}}^{T} = {D [(AA^{+})^{T}] AA^{+}}^{T} = [D(AA^{+})AA^{+}]^{T} ,

D(AA^{+}) = (1_{m} − AA^{+})D(A)A^{+} + [(1_{m} − AA^{+})D(A)A^{+}]^{T} ,

J_{2}(x_{2}) = −(1_{m} − AA^{+})D(A)A^{+}b − [D(A)A^{+}]^{T} (1_{m} − AA^{+})b + (1_{m} − AA^{+})D(b) =

= (1_{m} − AA^{+})J_{n−k}(x_{2}) − [D(A)A^{+}]^{T}f_{2}(x_{2}) .

¶

AA^{+}A = A , A^{+}AA^{+} = A^{+} ,

(AA^{+})^{T} = AA^{+} , (A^{+}A)^{T} = A^{+}A † quite generally, quadratic form

b^{T}AA^{+}b = b^{T}AA^{+}AA^{+}b = b^{T} (AA^{+})^{T} AA^{+}b = (A^{+}b)^{T} A^{T}AA^{+}b , A ∼ (m, n), b ∼ (m, 1) ∧ b ≠ 0, m ≥ n, is positive if and only if the matrix A^{T}A is positive definite and the vector A^{+}b is nonzero, i.e. if the matrix A has full rank

§ especially close to a local minimum the number of iterations is unbearably large

It is instantly evident that at least near the global minimum, where the l_{2}-norm of the vector f_{2}(x_{2}) is “small”, the importance of the second term decreases¶. It is just necessary to remember the commonly accepted approximation of the Hess matrix H(x) ∼ (n, n) of the function \frac{1}{2}∥f(x)∥^{2}

H(x) = \{\frac{∂^{2} \frac{1}{2}∥f(x)∥^{2}}{∂x_{i}∂x_{j}} \} , i = 1, …, n , j = 1, …, n ,

H(x) = J^{T}(x)J(x) + C(x)f(x) ,

C(x) = \{\frac{∂^{2}f_{i}}{∂x_{j}∂x_{k}} \} , i = 1, …, m , j = 1, …, n , k = 1, …, n ,

H(x) ≈ J^{T}(x)J(x) ,

which forms the basis of the Gauss-Newton method of solving the equation f(α) = 0. Hence, let us hope that the second (unpleasantly complicated) term −[D(A)A^{+}]^{T}f_{2}(x_{2}) in the relation for J_{2}(x_{2}) can be neglected without serious consequences, i.e. just the first term is enough to be taken into account^{2}

J_{2}(x_{2}) ≈ (1_{m} − AA^{+})J_{n−k}(x_{2})

and direct manipulation with the three-dimensional matrix D(A) of type (m, k, n − k) can be avoided (while both computer time and memory are saved). However, it should be noted that all the necessary derivatives have to be computed anyway, since they are “implicitly” included in the matrix J_{n−k}(x_{2}) = D(b) − D(A)A^{+}b†.

Thus, the described approximation is valuable especially for those, who want to \underline{easily and} \underline{quickly} modify almost any derivative algorithm of solving the equation f(α) = 0_{n} to a (possibly) more efficient algorithm of solving the “equivalent” equation f_{2}(α_{2}) = 0_{n−k} and are not willing to get over the loss of the already derived relations for calculation of the necessary partial derivatives. Since if we differentiate numerically, we only manage with n − k + 1 values (at different points) of the sum of squares of residuals ∥f_{2}(x_{2})∥^{2} in each iteration (to compute the gradient of the function \frac{1}{2}∥f_{2}(x_{2})∥^{2}) and n − k + 1 values of the function f_{2}(x_{2}) (to compute the Jacobi matrix J_{2}(x_{2})). The algorithm of calculation of ∥f_{2}(x_{2})∥^{2} becomes trivial after performing a suitable decomposition§ of the matrix A(x_{2}), the algorithm of calculation of f_{2}(x_{2}) is just a little more complicated. The question naturally arises if it is not more profitable to resign the use of analytical derivatives to the use of derivatives calculated according to common difference formulae.

If r(A) = k (i.e. if the matrix A(x_{2}) has full rank), then

A = Q^{T} ( \begin{matrix}U \\ 0 \end{matrix} ) = Q^{T} ( \begin{matrix}1_{k} \\ 0 \end{matrix} ) U = Q_{k}^{T}U ,
A^{+} = U^{−1}Q_{k} ,
f_{2}(x_{2}) = (1_{m} − Q_{k}^{T}Q_{k})b ,
∥f_{2}(x_{2})∥^{2} = ∥Q_{m−k}b∥^{2} ‡ ,

‡

f_{2}(x_{2}) = [Q^{T} ( \begin{matrix}1_{k} & 0 \\ 0 & 1_{m−k} \end{matrix} ) Q − Q^{T} ( \begin{matrix}1_{k} & 0 \\ 0 & 0 \end{matrix} ) Q] b = Q^{T} ( \begin{matrix}0 & 0 \\ 0 & 1_{m−k} \end{matrix} ) Qb = Q_{m−k}^{T}Q_{m−k}b ,

¶ However, the sum of squares of residuals can be large even at the global minimum if m > n.

† I will say once more, that the matrix J_{n−k}(x_{2}) is composed of the “right” submatrix of type (m, 1, n − k)

of the Jacobi matrix J(x) and vector x_{1} = A^{+}(x_{2})b(x_{2})

§ we only describe the QR- and singular value decomposition, which are the commonest

∥f_{2}(x_{2})∥^{2} = f_{2}^{T}(x_{2})f_{2}(x_{2}) = b^{T}Q_{m−k}^{T}Q_{m−k}Q_{m−k}^{T}Q_{m−k}b = b^{T}Q_{m−k}^{T}Q_{m−k}b = ∥Q_{m−k}b∥^{2}

J_{2}(x_{2}) = (1_{m} − Q_{k}^{T}Q_{k})J_{n−k}(x_{2}) − [D(A)U^{−1}Q_{k}]^{T}f_{2}(x_{2}) ,
J_{2}(x_{2}) ≈ (1_{m} − Q_{k}^{T}Q_{k})J_{n−k}(x_{2}) ,

Q ∼ (m, m) is an orthogonal matrix, U ∼ (k, k) is a regular upper-triangular matrix, or

A = Q^{T} ( \begin{matrix}U \\ 0 \end{matrix} ) P = Q^{T} ( \begin{matrix}1_{k} \\ 0 \end{matrix} ) UP = Q_{k}^{T}UP ,

A^{+} = P^{T}U^{−1}Q_{k} ,
J_{2}(x_{2}) = (1_{m} − Q_{k}^{T}Q_{k})J_{n−k}(x_{2}) − [D(A)P^{T}U^{−1}Q_{k}]^{T}f_{2}(x_{2}) ,

P ∼ (k, k) is a permutation matrix¶.

If r(A) ≡ r < k (i.e. if the matrix A(x_{2}) does not have full rank), then

A = Q^{T} ( \begin{matrix}U & V \\ 0 & 0 \end{matrix} ) P = Q^{T} ( \begin{matrix}1_{r} \\ 0 \end{matrix} ) ( U V ) P = Q_{r}^{T} ( U V ) P ,

A^{+} = P^{T} ( \begin{matrix}U^{−1} \\ 0 \end{matrix} ) Q_{r} = P^{T} ( \begin{matrix}1_{r} \\ 0 \end{matrix} ) U^{−1}Q_{r} = P_{r}^{T}U^{−1}Q_{r} ,

f_{2}(x_{2}) = (1_{m} − Q_{r}^{T}Q_{r})b ,
∥f_{2}(x_{2})∥^{2} = ∥Q_{m−r}b∥^{2} ,

J_{2}(x_{2}) = (1_{m} − Q_{r}^{T}Q_{r})J_{n−k}(x_{2}) − [D(A)P_{r}^{T}U^{−1}Q_{r}]^{T}f_{2}(x_{2}) ,

J_{2}(x_{2}) ≈ (1_{m} − Q_{r}^{T}Q_{r})J_{n−k}(x_{2}) ,

Q ∼ (m, m) is an orthogonal matrix, U ∼ (r, r) is a regular upper-triangular matrix, P ∼ (k, k) is a permutation matrix¶.

1. Golub G. H., Pereyra V., SIAM J. Numer. Anal. 10 (1973) 413

2. Maeder M., Zuberbühler A. D., Anal. Chem. 62 (1990) 2220

A^{+} = V^{T}Σ^{−1}U_{k} ,

A^{+} = V_{r}^{T}Σ^{−1}U_{r} ,

¶ Formally the same relations we get if we substitute the singular value decomposition for the QR-decompo-sition with column pivoting. If r(A) = k, then

A = U^{T} ( \begin{matrix}Σ \\ 0 \end{matrix} ) V = U^{T} ( \begin{matrix}1_{k} \\ 0 \end{matrix} ) ΣV = U_{k}^{T}ΣV , U ∼ (m, m) a V ∼ (k, k) are orthogonal matrices, Σ ≡ diag (σ_{1}, …, σ_{k}) ∼ (k, k) is a diagonal matrix with positive diagonal elements σ_{1} > … > σ_{k} > 0. If r(A) ≡ r < k, then

A = U^{T} ( \begin{matrix}Σ & 0 \\ 0 & 0 \end{matrix} ) V = U^{T} ( \begin{matrix}1_{r} \\ 0 \end{matrix} ) ( Σ 0 ) V = U_{r}^{T}Σ ( 1_{r} 0 ) V = U_{r}^{T}ΣV_{r} , U ∼ (m, m) a V ∼ (k, k) are orthogonal matrices, Σ ≡ diag (σ_{1}, …, σ_{r}) ∼ (r, r) is a diagonal matrix with positive diagonal elements σ_{1} > … > σ_{r} > 0.