<figure><img src="image_1.png" /></figure>

HAL
open science

# Solving shortest length least-squares problems via dynamic programming

Robert Kalaba, R Xu, W. Feng

► To cite this version:

Robert Kalaba, R Xu, W. Feng. Solving shortest length least-squares problems via dynamic programming. Journal of Optimization Theory and Applications, 1995, 85 (3), pp.613-632. 10.1007/BF02193059. hal-01406285

HAL Id: hal-01406285

https://hal.science/hal-01406285

Submitted on 1 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Solving Shortest Length Least-Squares Problems via Dynamic Programming<fnref n="1" />

R. KALABA,<fnref n="2" /> R. XU,<fnref n="3" /> AND W. FENG<fnref n="4" />

**Abstract.** If the matrix A is not of full rank, there may be many solutions to the problem of minimizing $||Ax - b||$ over x. Among such vectors x, the unique one for which $||x||$ is minimum is of importance in applications. This vector may be represented as $x = A^+b$. In this paper, the functional equation technique of dynamic programming is used to find the shortest solution to the least-squares problem in a sequential fashion. The algorithm is illustrated with an example.

**Key Words.** Dynamic programming, least-squares problems, principle of optimality, rank-deficient matrices.

## 1. Introduction

Consider the problem of finding the shortest length solution to the least-squares problem

$$Ax \approx b,$$

where A is an $m \times n$ matrix with rank r, the first r columns of which are independent, x is an $n \times 1$ vector, and b is an $m \times 1$ vector. Put in another way, we would like to find a vector x such that $x^T x$ is minimum among those that make Ax as close to b as possible in the sense that $(Ax - b)^T (Ax - b)$ is a minimum. As is well known, the solution to this problem is

$$x = A^+b,$$

<footnote n="2">²Professor, Departments of Biomedical Engineering, Electrical Engineering, and Economics, University of Southern California, Los Angeles, California.</footnote>

<footnote n="3">³PhD Candidate, School of Urban and Regional Planning, University of Southern California, Los Angeles, California.</footnote>

<footnote n="4">⁴PhD Student, Department of Economics, University of Southern California, Los Angeles, California.</footnote>

<footnote n="1">¹Our debt to the late Professor Richard Bellman is clear, and we wish to thank Professor Harriet Kagiwada for many stimulating conversations concerning least-squares problems over a long period of years.</footnote>
where $A^+$ is the Moore-Penrose generalized inverse of matrix A (Ref. 1). In this paper, we present an alternative approach to this shortest length least-squares problem (Ref. 2) via dynamic programming (Ref. 3).

Let $a_i$ be the ith column of the matrix A, and let $x_i$ be the ith element of the vector x, $i=1, 2, ..., n$. Then, the problem presented above is equivalent to finding the set of scalars, $x_1, x_2, ..., x_n$, such that $\|a_1x_1 + a_2x_2 + ... + a_nx_n - b\|$ is a minimum, and $x_1^2 + x_2^2 + ... + x_n^2$ is as small as possible.

In Section 2, an αQ-algorithm (Ref. 4) is introduced to solve the consistent set of linear algebraic equations

$$A_r x^{(r)} = c$$

in a sequential fashion; here, $A_r$ is the $m \times r$ submatrix of A whose columns are $a_1, a_2, ..., a_r$; $x^{(r)}$ is the $r \times 1$ subset of vector x whose elements are $x_1, x_2, ..., x_r$; and c is an $m \times 1$ vector. In Section 3, a βR-algorithm is introduced to find the set of values $x^{(r)}, x_{r+1}, x_{r+2}, ..., x_n$ such that $x_1^2 + x_2^2 + ... + x_n^2$ is minimum and the consistent set of linear equations

$$a_1 x_1 + a_2 x_2 + ... + a_n x_n = d$$

is satisfied, where d is an $m \times 1$ vector that is linearly dependent upon $a_1, a_2, ..., a_r$. Section 4 explains the connections between the above two problems and the original problem

$$Ax \approx b,$$

and an αQβR-algorithm is developed to solve the original problem. Section 5 provides an example to illustrate the procedures proposed. Discussions and extensions are given in Section 6, where we see that the αQβR-algorithm is capable of handling any shortest length least-squares problems as well as providing an algorithm for calculating $A^+$. Concluding remarks follow in Section 7.

## 2. The αQ-Algorithm

Consider the independent problem of finding a solution to the consistent set of linear algebraic equations

$$A_r x^{(r)} = c;$$

here, $A_r$ is an $m \times r$ matrix with rank $r$, whose columns are $a_1, a_2, ..., a_r$; $x^{(r)}$ is an $r \times 1$ vector whose elements are $x_1, x_2, ..., x_r$; and c is an $m \times 1$ vector. We use dynamic programming to determine the values of $x_1, x_2, ..., x_r$ sequentially.
Let $f_k(c_k)$ be the smallest square of the distance between $c_k$ and $\sum_{i=1}^k a_i x_i$, where $c_k$ is an $m \times 1$ vector and $k=1, 2, \dots, r$. Then, the Bellman principle of optimality leads to the recurrence relationship

$$f_k(c_k) = \min_{x_k} f_{k-1}(c_k - a_k x_k), \quad k=2, 3, \dots, r, \quad (1)$$

$$f_1(c_1) = \min_{x_1} (c_1 - a_1 x_1)^T (c_1 - a_1 x_1), \quad (2)$$

where, eventually,

$$c_r = c, \quad (3a)$$

$$c_k = c_{k+1} - a_{k+1} x_{k+1}^{\text{opt}}, \quad k=r-1, r-2, \dots, 1. \quad (3b)$$

Relation (2) is equivalent to

$$f_1(c_1) = \min_{x_1} (c_1^T c_1 - 2c_1^T a_1 x_1 + a_1^T a_1 x_1^2), \quad (4)$$

which is quadratic in $x_1$. Therefore, the first-order condition for the minimizing value of $x_1$ is

$$\frac{\partial}{\partial x_1} = 2a_1^T a_1 x_1 - 2a_1^T c_1 = 0. \quad (5)$$

It follows that

$$x_1^{\text{opt}} = a_1^T c_1 / a_1^T a_1, \quad a_1^T a_1 \neq 0, \quad (6)$$

or

$$x_1^{\text{opt}} = c_1^T a_1 / a_1^T a_1, \quad a_1^T a_1 \neq 0. \quad (7)$$

The inequalities in (6) and (7) hold, since $a_1$ is not the null vector. Substituting Eq. (6) or (7) into Eq. (4) gives

$$\begin{align*}
f_1(c_1) &= c_1^T c_1 - 2c_1^T (a_1 a_1^T / a_1^T a_1) c_1 + c_1^T (a_1 a_1^T / a_1^T a_1) c_1 \\
&= c_1^T c_1 - c_1^T (a_1 a_1^T / a_1^T a_1) c_1 \\
&= c_1^T (I - a_1 a_1^T / a_1^T a_1) c_1 \\
&= c_1^T Q_1 c_1,
\end{align*} \quad (8)$$

where

$$Q_1 = I - a_1 a_1^T / a_1^T a_1. \quad (9)$$

From Eqs. (8) and (9), we see that the minimum cost attainable when approximating $c_1$ by the first column of $A$, is $c_1^T Q_1 c_1$, where $Q_1$ is a symmetric $m \times m$ matrix.
Assuming that

$$f_{k-1}(\mathbf{c}_{k-1}) = \mathbf{c}_{k-1}^T \mathbf{Q}_{k-1} \mathbf{c}_{k-1},$$

where $\mathbf{Q}_{k-1}$ is a positive-definite symmetric $m \times m$ matrix, we now prove that $f_k(\mathbf{c}_k)$ has the form

$$f_k(\mathbf{c}_k) = \mathbf{c}_k^T \mathbf{Q}_k \mathbf{c}_k,$$

where $\mathbf{Q}_k$ is a symmetric $m \times m$ matrix and $k = 1, 2, \dots, r$.

From the recurrence relation (1), it can be seen that

$$ \begin{aligned} f_k(\mathbf{c}_k) &= \min_{x_k} f_{k-1}(\mathbf{c}_k - \mathbf{a}_k x_k) \\ &= \min_{x_k} \left[ (\mathbf{c}_k - \mathbf{a}_k x_k)^T \mathbf{Q}_{k-1} (\mathbf{c}_k - \mathbf{a}_k x_k) \right] \\ &= \min_{x_k} \left[ \mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k - 2 \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k x_k + \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k x_k^2 \right]. \end{aligned} \quad (10) $$

Realize that, in Eq. (10), the sum in the square brackets is quadratic in $x_k$. Therefore, the first-order condition for the minimizing value of $x_k$ is

$$\partial[\cdot]/\partial x_k = -2\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k + 2\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k x_k = 0. \quad (11)$$

Consequently,

$$x_k^{\text{opt}} = \frac{\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k}{\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k}, \quad \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k \neq 0, \quad (12)$$

or

$$x_k^{\text{opt}} = \frac{\mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k}{\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k}, \quad \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k \neq 0. \quad (13)$$

The last inequality will become clear later in Eqs. (41)–(44) and Eqs. (36)–(39), from which we can prove that $\mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k = \alpha_k^T \mathbf{a}_k \neq 0$. Substituting Eq. (12) or (13) into Eq. (10) gives

$$ \begin{aligned} f_k(\mathbf{c}_k) &= \mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k - 2\mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k x_k^{\text{opt}} + x_k^{\text{opt}} \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k x_k^{\text{opt}} \\ &= \mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k - 2\mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k (\mathbf{a}_k^T \mathbf{Q}_{k-1} / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k) \mathbf{c}_k \\ &\quad + \mathbf{c}_k^T (\mathbf{Q}_{k-1} \mathbf{a}_k / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k) \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k (\mathbf{a}_k^T \mathbf{Q}_{k-1} / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k) \mathbf{c}_k \\ &= \mathbf{c}_k^T \mathbf{Q}_{k-1} \mathbf{c}_k - \mathbf{c}_k^T (\mathbf{Q}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{Q}_{k-1} / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k) \mathbf{c}_k \\ &= \mathbf{c}_k^T (\mathbf{Q}_{k-1} - \mathbf{Q}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{Q}_{k-1} / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k) \mathbf{c}_k \\ &= \mathbf{c}_k^T \mathbf{Q}_k \mathbf{c}_k, \end{aligned} \quad (14) $$

where

$$\mathbf{Q}_k = \mathbf{Q}_{k-1} - \mathbf{Q}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{Q}_{k-1} / \mathbf{a}_k^T \mathbf{Q}_{k-1} \mathbf{a}_k, \quad k=2, 3, \dots, r. \quad (15)$$
Let

$$ \alpha_k = Q_{k-1} \alpha_k. \qquad (16) $$

Then, Eq. (15) is equivalent to

$$ Q_k = Q_{k-1} - \alpha_k \frac{\alpha_k^T}{\alpha_k^T \alpha_k}, \quad \alpha_k^T \alpha_k \neq 0, \quad k=2, 3, \dots, r. \qquad (17) $$

Define $\alpha_1 = \alpha_1$. Let us now examine the relationship between $\alpha_1, \alpha_2, \dots, \alpha_r$ and $\alpha_1, \alpha_2, \dots, \alpha_r$. It can be seen from Eq. (9) that, on the one hand, $Q_1$ premultiplying a vector that is a scalar multiple of $\alpha_1$ (say, $s_1\alpha_1$) yields

$$ \begin{aligned} Q_1 s_1 \alpha_1 &= (I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1}) s_1 \alpha_1 \\ &= s_1 \alpha_1 - s_1 (\alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} \alpha_1) = s_1 \alpha_1 - s_1 \alpha_1 = 0; \end{aligned} \qquad (18) $$

on the other hand, $Q_1$ premultiplying a vector (say, $v_1$) that is orthogonal to $\alpha_1$ yields

$$ Q_1 v_1 = (I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1}) v_1 = v_1 - 0 = v_1, \qquad (19) $$

for

$$ \alpha_1^T v_1 = 0. \qquad (20) $$

Thus, we conclude that any vector premultiplied by $Q_1$ equals the component of that vector that is perpendicular to $\alpha_1$. Consequently, according to the definition of $\alpha_2$ implied by Eq. (16), $\alpha_2$ is the component of $\alpha_2$ that is perpendicular to $\alpha_1$ or $\alpha_1$. That is,

$$ \alpha_2^T \alpha_1 = \alpha_1^T \alpha_2 = 0. \qquad (21) $$

Similarly, by Eqs. (17) and (9),

$$ Q_2 = I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} - \alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2}, \qquad (22) $$

it follows that, on the one hand, $Q_2$ premultiplying a vector that is a scalar multiple of $\alpha_1$ (say, $s_{21}\alpha_1$) or a scalar multiple of $\alpha_2$ (say, $s_{22}\alpha_2$) yields

$$ \begin{aligned} Q_2 s_{21} \alpha_1 &= (I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} - \alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2}) s_{21} \alpha_1 \\ &= s_{21} \alpha_1 - s_{21} (\alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} \alpha_1) - 0 \\ &= s_{21} \alpha_1 - s_{21} \alpha_1 = 0, \end{aligned} \qquad (23) $$

and

$$ \begin{aligned} Q_2 s_{22} \alpha_2 &= (I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} - \alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2}) s_{22} \alpha_2 \\ &= s_{22} \alpha_2 - 0 - s_{22} (\alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2} \alpha_2) \\ &= s_{21} \alpha_1 - s_{21} \alpha_1 = 0; \end{aligned} \qquad (24) $$
on the other hand, $\mathbb{Q}_2$ premultiplying a vector (say, $\mathbf{v}_2$) that is orthogonal to both $\mathbf{a}_1$ and $\mathbf{a}_2$ yields

$$ \mathbb{Q}_2 \mathbf{v}_2 = (\mathbf{I} - \mathbf{a}_1 \mathbf{a}_1^T / \mathbf{a}_1^T \mathbf{a}_1 - \mathbf{a}_2 \mathbf{a}_2^T / \mathbf{a}_2^T \mathbf{a}_2) \mathbf{v}_2 = \mathbf{v}_2 - 0 - 0 = \mathbf{v}_2, \quad (25) $$

for

$$ \mathbf{a}_1^T \mathbf{v}_2 = 0 \quad \text{and} \quad \mathbf{a}_2^T \mathbf{v}_2 = 0. \quad (26) $$

Thus, we conclude that any vector premultiplied by $\mathbb{Q}_2$ equals the component of that vector that is perpendicular to both $\mathbf{a}_1$ and $\mathbf{a}_2$. Consequently, based on the definition of $\mathbf{a}_3$ implied by Eq. (16), $\mathbf{a}_3$ is the component of $\mathbf{a}_3$ that is perpendicular to both $\mathbf{a}_1$ and $\mathbf{a}_2$. That is,

$$ \mathbf{a}_3^T \mathbf{a}_1 = \mathbf{a}_1^T \mathbf{a}_3 = 0 \quad \text{and} \quad \mathbf{a}_3^T \mathbf{a}_2 = \mathbf{a}_2^T \mathbf{a}_3 = 0. \quad (27) $$

This way, it can be inferred that $\mathbf{a}_k$ is actually the component of $\mathbf{a}_k$ that is orthogonal to all the earlier $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_{k-1}$, and $k=1, 2, \dots, r$. Put in another way,

$$ \mathbf{a}_k^T \mathbf{a}_l = 0, \quad \text{for all } k \neq l. \quad (28) $$

Furthermore, $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ and $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ span exactly the same vector space. This is so because $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ can be expressed by a linear combination of $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ and vice versa. That is,

$$ \mathbf{a}_1 = \mathbf{a}_1, \quad (29) $$

$$ \mathbf{a}_1 = \mathbf{a}_1, \quad (30) $$

$$ \mathbf{a}_2 = \mathbf{a}_2 + u_{12} \mathbf{a}_1, \quad (31) $$

$$ \mathbf{a}_2 = \mathbf{a}_2 - u_{12} \mathbf{a}_1 = \mathbf{a}_2 - u_{12} \mathbf{a}_1, \quad (32) $$

$$ \mathbf{a}_3 = \mathbf{a}_3 + u_{23} \mathbf{a}_2 + u_{13} \mathbf{a}_1, \quad (33) $$

$$ \begin{aligned} \mathbf{a}_3 &= \mathbf{a}_3 - u_{23} \mathbf{a}_2 - u_{13} \mathbf{a}_1 = \mathbf{a}_3 - u_{23} (\mathbf{a}_2 - u_{12} \mathbf{a}_1) - u_{13} \mathbf{a}_1 \\ &= \mathbf{a}_3 - u_{23} \mathbf{a}_2 + (u_{23} u_{12} - u_{13}) \mathbf{a}_1, \quad (34) \end{aligned} $$

$\vdots$

where the $u$'s are all scalars. In matrix form, we have

$$ [\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r] = [\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r] \begin{bmatrix} 1 & u_{12} & u_{13} & \dots & u_{1r} \\ 0 & 1 & u_{23} & \dots & u_{2r} \\ 0 & 0 & 1 & \dots & u_{3r} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \end{bmatrix}. \quad (35) $$
It follows that

$$ \alpha_1^T \alpha_1 = \alpha_1^T \alpha_1, \qquad (36) $$

$$ \alpha_2^T \alpha_2 = \alpha_2^T (\alpha_2 + u_{12}\alpha_1) = \alpha_2^T \alpha_2 + 0 = \alpha_2^T \alpha_2, \qquad (37) $$

$$ \alpha_3^T \alpha_3 = \alpha_3^T (\alpha_3 + u_{23}\alpha_2 + u_{13}\alpha_1) = \alpha_3^T \alpha_3 + 0 + 0 = \alpha_3^T \alpha_3, \qquad (38) $$

∴

$$ \alpha_r^T \alpha_r = \alpha_r^T (\alpha_r + u_{(r-1)r} \alpha_{r-1} + \dots + u_{2r} \alpha_2 + u_{1r} \alpha_1) = \alpha_r^T \alpha_r. \qquad (39) $$

Substituting Eqs. (36)-(39) into Eq. (17) provides

$$ Q_k = Q_{k-1} - \alpha_k \frac{\alpha_k^T}{\alpha_k^T \alpha_k}, \quad \alpha_k^T \alpha_k \neq 0, \quad k=2, 3, \dots, r. \qquad (40) $$

Define $Q_0 = I$. Then, Eq. (40) holds for $k=1, 2, \dots, r$. More specifically,

$$ Q_0 = I, \qquad (41) $$

$$ Q_1 = I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1}, \quad \alpha_1^T \alpha_1 \neq 0, \qquad (42) $$

$$ Q_2 = I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} - \alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2}, \qquad (43a) $$

$$ \alpha_1^T \alpha_1 \neq 0 \quad \text{and} \quad \alpha_2^T \alpha_2 \neq 0, \qquad (43b) $$

∴

$$ Q_r = I - \alpha_1 \frac{\alpha_1^T}{\alpha_1^T \alpha_1} - \alpha_2 \frac{\alpha_2^T}{\alpha_2^T \alpha_2} - \dots - \alpha_r \frac{\alpha_r^T}{\alpha_r^T \alpha_r}, \qquad (44a) $$

$$ \alpha_1^T \alpha_1 \neq 0, \quad \alpha_2^T \alpha_2 \neq 0, \quad \dots, \quad \alpha_r^T \alpha_r \neq 0. \qquad (44b) $$

The reason $\alpha_k^T \alpha_k \neq 0$ is because of the assumption that the columns of $A_r$ are independent; hence, each of them does have a component that is perpendicular to all its proceeding column vectors. This is true for all $k, k=1, 2, \dots, r$. Substituting Eqs. (36)-(39) into Eqs. (6) and (12) provides

$$ x_1^{\text{opt}} = \frac{\alpha_1^T c_1}{\alpha_1^T \alpha_1}, \quad \alpha_1^T \alpha_1 \neq 0, \qquad (45) $$

and

$$ x_k^{\text{opt}} = \frac{\alpha_k^T c_k}{\alpha_k^T \alpha_k}, \quad \alpha_k^T \alpha_k \neq 0. \qquad (46) $$

These are the simplified formulas for the optimal set of scalars $x_1, x_2, \dots, x_r$.

Notice that $\alpha_1, \alpha_2, \dots, \alpha_r$ are independent; hence, none of them is the null vector. By assumption, the system

$$ A_r x^{(r)} = c $$
is consistent. It follows that there is one and only one set of scalars $x_1, x_2, \dots, x_r$ that satisfy

$$ \mathbf{a}_1\mathbf{x}_1 + \mathbf{a}_2\mathbf{x}_2 + \dots + \mathbf{a}_r\mathbf{x}_r = \mathbf{c}. $$

Clearly, this set of scalars $x_1, x_2, \dots, x_r$ is the shortest available one.

In short, to find the solution to the consistent system

$$ \mathbf{A}, \mathbf{x}^{(r)} = \mathbf{c} $$

sequentially, we may use the following $\alpha Q$-algorithm. We first carry out the sequence of calculations

$$ \mathbf{a}_1 = \mathbf{a}_1, \qquad (47a) $$

$$ \mathbf{Q}_1 = I - \mathbf{a}_1 \mathbf{a}_1^T / \mathbf{a}_1^T \mathbf{a}_1, \quad \mathbf{a}_1^T \mathbf{a}_1 \neq 0, \qquad (47b) $$

$$ \mathbf{a}_2 = \mathbf{Q}_1 \mathbf{a}_2, \qquad (48a) $$

$$ \mathbf{Q}_2 = \mathbf{Q}_1 - \mathbf{a}_2 \mathbf{a}_2^T / \mathbf{a}_2^T \mathbf{a}_2, \quad \mathbf{a}_2^T \mathbf{a}_2 \neq 0, \qquad (48b) $$

$$ \mathbf{a}_r = \mathbf{Q}_{r-1} \mathbf{a}_r, \qquad (49a) $$

$$ \mathbf{Q}_r = \mathbf{Q}_{r-1} - \mathbf{a}_r \mathbf{a}_r^T / \mathbf{a}_r^T \mathbf{a}_r, \quad \mathbf{a}_r^T \mathbf{a}_r \neq 0. \qquad (49b) $$

The set of scalars $x_1, x_2, \dots, x_r$ that satisfy

$$ \mathbf{A}_r \mathbf{x}^{(r)} = \mathbf{c} $$

is then obtained by

$$ x_r^{\text{opt}} = \mathbf{a}_r^T \mathbf{c} / \mathbf{a}_r^T \mathbf{a}_r, \quad \mathbf{a}_r^T \mathbf{a}_r \neq 0, \qquad (50) $$

$$ x_{r-1}^{\text{opt}} = \mathbf{a}_{r-1}^T (\mathbf{c} - \mathbf{a}_r x_r^{\text{opt}}) / \mathbf{a}_{r-1}^T \mathbf{a}_{r-1}, \quad \mathbf{a}_{r-1}^T \mathbf{a}_{r-1} \neq 0, \qquad (51) $$

$$ x_1^{\text{opt}} = \mathbf{a}_1^T \left( \mathbf{c} - \sum_{i=2}^{r} \mathbf{a}_i x_i^{\text{opt}} \right) / \mathbf{a}_1^T \mathbf{a}_1, \quad \mathbf{a}_1^T \mathbf{a}_1 \neq 0. \qquad (52) $$

### 3. The βR-Algorithm

Consider the independent problem of finding the shortest length solution to the consistent set of linear algebraic equations

$$ \mathbf{A}\mathbf{x} = \mathbf{d}; $$

here, $\mathbf{A}$ is an $m \times n$ matrix of rank $r$ having its first $r$ columns $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ linearly independent; $\mathbf{x}$ is an $n \times 1$ vector; and $\mathbf{d}$ is an $m \times 1$ vector. Put in
another way, we want to find the set of scalars $x_1, x_2, \dots, x_n$ such that the assumed consistent set of linear equations

$$ \mathbf{a}_1 x_1 + \mathbf{a}_2 x_2 + \dots + \mathbf{a}_n x_n = \mathbf{d} $$

is satisfied and $x_1^2 + x_2^2 + \dots + x_n^2$ is minimum. Again, we use dynamic programming to solve this problem sequentially.

Let $g_k(\mathbf{d}_k)$ be the value of $x_1^2 + x_2^2 + \dots + x_k^2$ when using an optimal set of scalars $x_1, x_2, \dots, x_k$; here,

$$ \mathbf{a}_1 x_1 + \mathbf{a}_2 x_2 + \dots + \mathbf{a}_k x_k = \mathbf{d}_k $$

is a consistent set of linear algebraic equations; $\mathbf{d}_k$ is an $m \times 1$ vector; and $k=r, r+1, \dots, n$. Optimal is in the sense that this particular set of scalars $x_1, x_2, \dots, x_k$ renders the value of $x_1^2 + x_2^2 + \dots + x_k^2$ a minimum. Then, the Bellman principle of optimality leads to the recurrence relationship

$$ g_k(\mathbf{d}_k) = \min_{x_k} \{x_k^2 + g_{k-1}(\mathbf{d}_k - \mathbf{a}_k x_k)\}, \quad k=r+1, r+2, \dots, n, \quad (53) $$

where, eventually,

$$ \mathbf{d}_k = \mathbf{d}_{k+1} - \mathbf{a}_{k+1} x_{k+1}, \quad k=n-1, n-2, \dots, r, \quad (54a) $$

$$ \mathbf{d}_n = \mathbf{d}. \quad (54b) $$

In Eq. (53), the reason for which $k$ starts from $r+1$ is that there is one and only one set of scalars $x_1, x_2, \dots, x_r$ that fulfills the constraint

$$ \mathbf{a}_1 x_1 + \mathbf{a}_2 x_2 + \dots + \mathbf{a}_r x_r = \mathbf{d}_r; $$

hence, there is no freedom in choosing the scalars $x_1, x_2, \dots, x_k$ when $k \le r$.

Let $\mathbf{A}_r$ be the $m \times r$ matrix whose columns are $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$, and let $\mathbf{x}^{(r)}$ be the $r \times 1$ vector whose elements are $x_1, x_2, \dots, x_r$. Since $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ are independent, the matrix $\mathbf{A}_r^T \mathbf{A}_r$ is of dimension $r \times r$ with rank $r$; hence, its inverse $(\mathbf{A}_r^T \mathbf{A}_r)^{-1}$ exists. Thus, the unique solution to the consistent set of linear algebraic equations

$$ \mathbf{A}_r \mathbf{x}^{(r)} = \mathbf{d}_r, \quad (54c) $$

in vector form, is

$$ \mathbf{x}^{(r)} = (\mathbf{A}_r^T \mathbf{A}_r)^{-1} \mathbf{A}_r^T \mathbf{d}_r. \quad (55) $$

Consequently, by definition,

$$ \begin{align} g_r(\mathbf{d}_r) &= x_1^2 + x_2^2 + \dots + x_r^2 = \mathbf{x}^{(r)T} \mathbf{x}^{(r)} \nonumber \\ &= \mathbf{d}_r^T \mathbf{A}_r [ (\mathbf{A}_r^T \mathbf{A}_r)^{-1} ]^T (\mathbf{A}_r^T \mathbf{A}_r)^{-1} \mathbf{A}_r^T \mathbf{d}_r \nonumber \\ &= \mathbf{d}_r^T \mathbf{A}_r (\mathbf{A}_r^T \mathbf{A}_r)^{-1} (\mathbf{A}_r^T \mathbf{A}_r)^{-1} \mathbf{A}_r^T \mathbf{d}_r. \tag{56} \end{align} $$
Denote

$$ \mathbf{R}_r = \mathbf{A}_r (\mathbf{A}_r^T \mathbf{A}_r)^{-1} (\mathbf{A}_r^T \mathbf{A}_r)^{-1} \mathbf{A}_r^T. \quad (57) $$

Then,

$$ g_r(\mathbf{d}_r) = \mathbf{d}_r^T \mathbf{R}_r \mathbf{d}_r, \quad (58) $$

where $\mathbf{R}_r$ is an $m \times m$ symmetric matrix. Later, in Section 6, we will give an alternative formula for $\mathbf{R}_r$ with no matrix inversion involved.

Assuming that

$$ g_{k-1}(\mathbf{d}_{k-1}) = \mathbf{d}_{k-1}^T \mathbf{R}_{k-1} \mathbf{d}_{k-1}, $$

where $\mathbf{R}_{k-1}$ is an $m \times m$ positive-definite symmetric matrix, we now prove that $g_k(\mathbf{d}_k)$ has the form

$$ g_k(\mathbf{d}_k) = \mathbf{d}_k^T \mathbf{R}_k \mathbf{d}_k, $$

where $\mathbf{R}_k$ is an $m \times m$ symmetric matrix and $k=r+1, r+2, \dots, n$.

From the recurrence relation (53), we see that, for $k=r+1, r+2, \dots, n$,

$$ \begin{aligned} g_k(\mathbf{d}_k) &= \min_{x_k} \{x_k^2 + g_{k-1}(\mathbf{d}_k - \mathbf{a}_k x_k)\} \\ &= \min_{x_k} \{x_k^2 + (\mathbf{d}_k - \mathbf{a}_k x_k)^T \mathbf{R}_{k-1} (\mathbf{d}_k - \mathbf{a}_k x_k)\} \\ &= \min_{x_k} \{x_k^2 + \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{d}_k - 2 \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{a}_k x_k + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k x_k^2\} \\ &= \min_{x_k} \{(1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k) x_k^2 + \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{d}_k - 2 \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{a}_k x_k\}. \end{aligned} \quad (59) $$

Since the first-order condition for the minimizing value of $x_k$ is

$$ \partial\{\cdot\} / \partial x_k = 2(1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k) x_k - 2 \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{a}_k = 0, \quad (60) $$

we find that

$$ x_k^{\text{opt}} = \frac{\mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{a}_k}{1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k}, \quad (61) $$

or

$$ x_k^{\text{opt}} = \frac{\mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{d}_k}{1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k}. \quad (62) $$

Substituting Eq. (61) or Eq. (62) into Eq. (59) gives

$$ \begin{aligned} g_k(\mathbf{d}_k) &= (1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k) [x_k^{\text{opt}}]^2 - 2 \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{d}_k x_k^{\text{opt}} + \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{d}_k \\ &= (1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k) [\mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{d}_k / (1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k)^2] \\ &\quad - 2 \mathbf{d}_k^T [\mathbf{R}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{R}_{k-1} / (1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k)] \mathbf{d}_k + \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{d}_k \\ &= \mathbf{d}_k^T \mathbf{R}_{k-1} \mathbf{d}_k - \mathbf{d}_k^T [\mathbf{R}_{k-1} \mathbf{a}_k \mathbf{a}_k^T \mathbf{R}_{k-1} / (1 + \mathbf{a}_k^T \mathbf{R}_{k-1} \mathbf{a}_k)] \mathbf{d}_k \end{aligned} $$
$$ 
\begin{aligned}
&= d_k^T [R_{k-1} - R_{k-1} a_k a_k^T R_{k-1} / (1 + a_k^T R_{k-1} a_k)] d_k \\
&= d_k^T R_k d_k,
\end{aligned} 
\quad (63) 
$$

where

$$ R_k = R_{k-1} - \frac{R_{k-1} a_k a_k^T R_{k-1}}{1 + a_k^T R_{k-1} a_k}, \quad k=r+1, r+2, \dots, n. \quad (64) $$

Notice that the denominator,

$$ 1 + a_k^T R_{k-1} a_k \neq 0, $$

always holds. Let

$$ \beta_k = R_{k-1} a_k, \quad k=r+1, r+2, \dots, n. \quad (65) $$

Equation (64) is equivalent to

$$ R_k = R_{k-1} - \frac{\beta_k \beta_k^T}{1 + a_k^T \beta_k}, \quad k=r+1, r+2, \dots, n. \quad (66) $$

Substituting Eq. (66) into Eq. (61) or Eq. (62) gives

$$ x_k^{\text{opt}} = \frac{d_k^T \beta_k}{1 + a_k^T \beta_k}, \quad (67) $$

or

$$ x_k^{\text{opt}} = \frac{\beta_k^T d_k}{1 + a_k^T \beta_k}, \quad (68) $$

where

$$ d_n = d, \quad (69a) $$

$$ d_k = d_{k+1} - a_{k+1} x_{k+1}^{\text{opt}}, \quad k=n-1, n-2, \dots, r. \quad (69b) $$

Therefore, the sequence of $R_r$, $R_{r+1}$, ..., $R_n$ is now obtainable through formulas (57) and (66). The optimal values of $x^{(r)}$, $x_{r+1}$, $x_{r+2}$, ..., $x_n$ can be calculated via Eq. (55) and Eq. (67) or Eq. (68).

In brief, to find the shortest length solution to the consistent set of linear algebraic equations

$$ Ax = d, $$

we use the following $\beta R$-algorithm. We carry out the sequence of calculations

$$ R_r = A_r (A_r^T A_r)^{-1} (A_r^T A_r)^{-1} A_r^T, \quad (70) $$

$$ \beta_{r+1} = R_r a_{r+1}, \quad (71a) $$

$$ R_{r+1} = R_r - \frac{\beta_{r+1} R_{r+1}^T}{1 + a_{r+1}^T \beta_{r+1}}, \quad (71b) $$

$$ \beta_{r+2} = R_{r+1} a_{r+2}, \quad (72a) $$
$$ \mathbf{R}_{r+2} = \mathbf{R}_{r+1} - \frac{\beta_{r+2}\beta_{+2}^T}{1 + \mathbf{a}_{r+2}^T\beta_{r+2}}, \quad (72b) $$

$$ \beta_n = \mathbf{R}_{n-1}\mathbf{a}_n, \quad (73a) $$

$$ \mathbf{R}_n = \mathbf{R}_{n-1} - \frac{\beta_n \beta_n^T}{1 + \mathbf{a}_n^T \beta_n}. \quad (73b) $$

The optimal set of scalars $x_1, x_2, \dots, x_n$ is then obtained by

$$ x_n^{\text{opt}} = \frac{\beta_n^T \mathbf{d}}{1 + \mathbf{a}_n^T \beta_n}, \quad (74) $$

$$ x_{n-1}^{\text{opt}} = \frac{\beta_{n-1}^T (\mathbf{d} - \mathbf{a}_n x_n^{\text{opt}})}{1 + \mathbf{a}_{n-1}^T \beta_{n-1}}, \quad (75) $$

$$ x_{r+1}^{\text{opt}} = \frac{\beta_{r+1}^T \left( \mathbf{d} - \sum_{i=r+2}^{n} \mathbf{a}_i x_i^{\text{opt}} \right)}{1 + \mathbf{a}_{r+1}^T \beta_{r+1}}, \quad (76) $$

$$ \begin{bmatrix} x_1^{\text{opt}} \\ x_2^{\text{opt}} \\ \vdots \\ x_r^{\text{opt}} \end{bmatrix} = (A_r^T A_r)^{-1} A_r^T \left( \mathbf{d} - \sum_{i=r+1}^{n} \mathbf{a}_i x_i^{\text{opt}} \right). \quad (77) $$

### 4. The αQβR-Procedure

We now turn to the original least-squares problem of finding the shortest length solution to

$$ \mathbf{A}\mathbf{x} \approx \mathbf{b}; $$

here, A is an $m \times n$ matrix of rank $r$, having its first $r$ columns, $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$, linearly independent; x is an $n \times 1$ vector whose elements are $x_1, x_2, \dots, x_n$; and b is an $m \times 1$ vector. The αQ and βR-algorithms can be extended and used here.

Recalling that, in carrying out the αQ-algorithm, we have also obtained $\mathbf{a}_k$, the component of the $k$th column of the matrix A, that is perpendicular to all its preceeding column vectors $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_{k-1}$ and $k=1, 2, \dots, r$. It can be inferred that the component of b that is perpendicular to all the vectors $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$ can be obtained by

$$ \mathbf{a}_{r+1} = \mathbf{Q}_r \mathbf{b}. \quad (78) $$

Since A is the augmented matrix of A, discussed in Section 3, and since each of the additional columns $\mathbf{a}_{r+1}, \mathbf{a}_{r+2}, \dots, \mathbf{a}_n$, in A, can be expressed by a
certain linear combination of $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_r$, it follows that $\mathbf{a}_{r+1}$ is also the component of $\mathbf{b}$ that is perpendicular to all the columns $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n$. Or alternatively, $\mathbf{b} - \mathbf{a}_{r+1}$ is the component of $\mathbf{b}$ that can be expressed by a linear combination of $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n$. Therefore,

$$ \mathbf{a}_1 x_1 + \mathbf{a}_2 x_2 + \dots + \mathbf{a}_n x_n = \mathbf{b} - \mathbf{a}_{r+1} $$

is a consistent set of linear algebraic equations. The problem of finding the shortest solution to the consistent system of equations

$$ \mathbf{A}\mathbf{x} = \mathbf{b} - \mathbf{a}_{r+1} $$

is in the form that can be solved by the βR-algorithm introduced in Section 4.

In summary, the general procedure (αQβR-algorithm) proposed to find the shortest solution to the least-squares problem $\mathbf{A}\mathbf{x} \approx \mathbf{b}$ is as follows:

(i) Calculate $\mathbf{a}_1, \mathbf{Q}_1, \mathbf{a}_2, \mathbf{Q}_2, \dots, \mathbf{a}_r, \mathbf{Q}_r$ by Eqs. (47)-(49).

(ii) Calculate $\mathbf{a}_{r+1}$ by Eq. (78).

(iii) Calculate $\mathbf{R}_r, \mathbf{\beta}_{r+1}, \mathbf{R}_{r+1}, \dots, \mathbf{\beta}_{n-1}, \mathbf{R}_{n-1}, \mathbf{\beta}_n$ by Eqs. (70)-(73).

(iv) Calculate the optimal set of scalars $x_1, x_2, \dots, x_n$ by

$$ \mathbf{b}_n = \mathbf{b} - \mathbf{a}_{r+1}, \qquad (79a) $$

$$ x_n^{\text{opt}} = \frac{\mathbf{\beta}_n^T \mathbf{b}_n}{1 + \mathbf{a}_n^T \mathbf{\beta}_n}, \qquad (79b) $$

$$ \mathbf{b}_{n-1} = \mathbf{b}_n - \mathbf{a}_n x_n^{\text{opt}}, \qquad (80a) $$

$$ x_{n-1}^{\text{opt}} = \frac{\mathbf{\beta}_{n-1}^T \mathbf{b}_{n-1}}{1 + \mathbf{a}_{n-1}^T \mathbf{\beta}_{n-1}}, \qquad (80b) $$

$$ \mathbf{b}_{r+1} = \mathbf{b}_{r+2} - \mathbf{a}_{r+2} x_{r+2}^{\text{opt}}, \qquad (81a) $$

$$ x_{r+1}^{\text{opt}} = \frac{\mathbf{\beta}_{r+1}^T \mathbf{b}_{r+1}}{1 + \mathbf{a}_{r+1}^T \mathbf{\beta}_{r+1}}, \qquad (81b) $$

$$ \mathbf{b}_r = \mathbf{b}_{r+1} - \mathbf{a}_{r+1} x_{r+1}^{\text{opt}}, \qquad (82a) $$

$$ x_r^{\text{opt}} = \frac{\mathbf{a}_r^T \mathbf{b}_r}{\mathbf{a}_r^T \mathbf{a}_r}, \quad \mathbf{a}_r^T \mathbf{a}_r \neq 0, \qquad (82b) $$

$$ \mathbf{b}_{r-1} = \mathbf{b}_r - \mathbf{a}_r x_r^{\text{opt}}, \qquad (83a) $$
$$x_{r-1}^{\text{opt}} = \frac{\mathbf{a}_{r-1}^{T} \mathbf{b}_{r-1}}{\mathbf{a}_{r-1}^{T} \mathbf{a}_{r-1}}, \quad \mathbf{a}_{r-1}^{T} \mathbf{a}_{r-1} \neq 0, \qquad (83b)$$

$$\mathbf{b}_1 = \mathbf{b}_2 - \mathbf{a}_2 x_2^{\text{opt}}, \qquad (84a)$$

$$x_1^{\text{opt}} = \frac{\mathbf{a}_1^T \mathbf{b}_1}{\mathbf{a}_1^T \mathbf{a}_1}, \quad \mathbf{a}_1^T \mathbf{a}_1 \neq 0. \qquad (84b)$$

## 5. Example

Let us consider an example, in which we want to find the shortest length vector x,

$$\mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix},$$

such that $||\mathbf{Ax} - \mathbf{b}||$ is minimum, where

$$\mathbf{A} = \begin{bmatrix} 1 & 1 & 2 & 2 \\ 0 & 1 & 1 & 2 \\ 0 & 1 & 1 & 2 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} 2 \\ -1 \\ 1 \end{bmatrix}.$$

Apparently, A has rank $r=2$. Its first two columns are independent, and its last two columns are certain linear combinations of the first two columns. Denote the columns of A by vectors

$$\mathbf{a}_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{a}_2 = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}, \quad \mathbf{a}_3 = \begin{bmatrix} 2 \\ 1 \\ 1 \end{bmatrix}, \quad \mathbf{a}_4 = \begin{bmatrix} 2 \\ 2 \\ 2 \end{bmatrix}.$$

Let

$$\mathbf{A}_2 = \begin{bmatrix} 1 & 1 \\ 0 & 1 \\ 0 & 1 \end{bmatrix}.$$

We follow the procedures prescribed in Section 4 to solve this shortest length least-squares problem.
(i) Calculate $\alpha_1$, $Q_1$, $\alpha_2$, $Q_2$ by Eqs. (47)–(49). We have

$$ \alpha_1 = \alpha_1 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \qquad (85) $$

$$ Q_1 = I - \alpha_1 \alpha_1^T / \alpha_1^T \alpha_1 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \qquad (86) $$

$$ \alpha_2 = Q_1 \alpha_2 = \begin{bmatrix} 0 \\ 1 \\ 1 \end{bmatrix}, \qquad (87) $$

$$ Q_2 = Q_1 - \alpha_2 \alpha_2^T / \alpha_2^T \alpha_2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1/2 & -1/2 \\ 0 & -1/2 & 1/2 \end{bmatrix}. \qquad (88) $$

(ii) Calculate $\alpha_3$ by Eq. (78). We have

$$ \alpha_3 = Q_2 b = \begin{bmatrix} 0 \\ -1 \\ 1 \end{bmatrix}. \qquad (89) $$

(iii) Calculate $R_2$, $\beta_3$, $R_3$, $\beta_4$ by Eqs. (70)–(73). Since

$$ A_2^T A_2 = \begin{bmatrix} 1 & 1 \\ 1 & 3 \end{bmatrix}, \qquad (90) $$

it follows that

$$ (A_2^T A_2)^{-1} = \begin{bmatrix} 3/2 & -1/2 \\ -1/2 & 1/2 \end{bmatrix}. \qquad (91) $$

Therefore,

$$ R_2 = A_2 (A_2^T A_2)^{-1} (A_2^T A_2)^{-1} A_2^T = \begin{bmatrix} 1 & -1/2 & -1/2 \\ -1/2 & 1/2 & 1/2 \\ -1/2 & 1/2 & 1/2 \end{bmatrix}, \qquad (92) $$

$$ \beta_3 = R_2 \alpha_3 = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \qquad (93) $$
$$ \mathbf{R}_3 = \mathbf{R}_2 - \frac{\beta_3 \beta_3^T}{1 + \mathbf{a}_3^T \beta_3} = \begin{bmatrix} \frac{2}{3} & -\frac{1}{2} & -\frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{bmatrix}, \quad (94) $$

$$ \beta_4 = \mathbf{R}_3 \mathbf{a}_4 = \begin{bmatrix} -2/3 \\ 1 \\ 1 \end{bmatrix}. \quad (95) $$

(iv) Calculate the optimal set of scalars $x_4, x_3, x_2, x_1$ by Eqs. (79)–(84). We have

$$ \mathbf{b}_4 = \mathbf{b} - \alpha_3 = \begin{bmatrix} 2 \\ 0 \\ 0 \end{bmatrix}, \quad (96) $$

$$ x_4^{\text{opt}} = \beta_4^T \mathbf{b}_4 / (1 + \alpha_4^T \beta_4) = -4/11, \quad (97) $$

$$ \mathbf{b}_3 = \mathbf{b}_4 - x_4^{\text{opt}} \mathbf{a}_4 = \begin{bmatrix} 30/11 \\ 8/11 \\ 8/11 \end{bmatrix}, \quad (98) $$

$$ x_3^{\text{opt}} = \frac{\beta_3^T \mathbf{b}_3}{1 + \alpha_3^T \beta_3} = 10/11, \quad (99) $$

$$ \mathbf{b}_2 = \mathbf{b}_3 - x_3^{\text{opt}} \mathbf{a}_3 = \begin{bmatrix} 10/11 \\ -2/11 \\ -2/11 \end{bmatrix}, \quad (100) $$

$$ x_2^{\text{opt}} = \alpha_2^T \mathbf{b}_2 / \alpha_2^T \alpha_2 = -2/11, \quad (101) $$

$$ \mathbf{b}_1 = \mathbf{b}_2 - x_2^{\text{opt}} \mathbf{a}_2 = \begin{bmatrix} 12/11 \\ 0 \\ 0 \end{bmatrix}, \quad (102) $$

$$ x_1^{\text{opt}} = \frac{\alpha_1^T \mathbf{b}_1}{\alpha_1^T \alpha_1} = \frac{12}{11}. \quad (103) $$

Therefore, the solution to this problem is

$$ \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} \frac{12}{11} \\ -\frac{2}{11} \\ \frac{10}{11} \\ -\frac{4}{11} \end{bmatrix}, $$

which can be verified to be correct by the elimination method or otherwise.
## 6. Discussion and Extension

The procedure prescribed in Section 4 can also be modified to handle a more general shortest length least-squares problem

$$Ax \cong b,$$

where A is an $m \times n$ matrix whose rank is unknown. The modifications include finding the rank and rearranging the columns of A so that it will have the form described in Section 1. These can be done via the $\alpha Q$-algorithm. More specifically, we calculate the values of $\alpha_i$ and $Q_i$ for all the columns of A, namely, for $i=1, 2, ..., n$. If for any column $l$, we encounter $\alpha_l=0$, we will move this column to the right end of A and continue. This way, after carrying out the $\alpha Q$-algorithm, all the columns with $\alpha_l \neq 0$ are at the left, and all the columns with $\alpha_l=0$ are at the right. The number of the columns with $\alpha_l \neq 0$ is then the rank of A.

It is also worthwhile to point out that there is an alternative way to calculate the $m \times m$ matrix $R_r$ in Eq. (57). Recall that, in Eq. (35), the $m \times r$ matrix $A_r$ is factored as

$$A_r = \alpha U, \qquad (104)$$

where

$$\alpha = [\alpha_1, \alpha_2, ..., \alpha_r]$$

and U is an $r \times r$ nonsingular upper triangular matrix with 1's on the main diagonal. Since

$$\alpha_2 = \alpha_2 + u_{12}\alpha_1, \qquad (105)$$

it follows that premultiplying both sides of Eq. (105) by $\alpha_i^T$ yields

$$\alpha_1^T \alpha_2 = \alpha_1^T \alpha_2 + u_{12} \alpha_1^T \alpha_1. \qquad (106)$$

Furthermore, since $\alpha_1^T \alpha_2 = 0$, Eq. (105) is equivalent to

$$\alpha_1^T \alpha_2 = u_{12} \alpha_1^T \alpha_1, \qquad (107)$$

or

$$u_{12} = \alpha_1^T \alpha_2 / \alpha_1^T \alpha_1. \qquad (108)$$

Similarly, since

$$\alpha_3 = \alpha_3 + u_{23}\alpha_2 + u_{13}\alpha_1, \qquad (109)$$
premultiplying both sides of Eq. (109) by $\alpha_2^T$ and $\alpha_1^T$ produces

$$u_{23} = \alpha_2^T \alpha_3 / \alpha_2^T \alpha_2 \qquad (110)$$

and

$$u_{13} = \alpha_1^T \alpha_3 / \alpha_1^T \alpha_1, \qquad (111)$$

respectively. This way, it can be shown that the elements $u_{ij}$, in the matrix U, can be obtained by the following formulas:

$$u_{ij} = \alpha_i^T \alpha_j / \alpha_i^T \alpha_i, \quad i < j, \qquad (112)$$

$$u_{ij} = 1, \quad i = j, \qquad (113)$$

$$u_{ij} = 0, \quad i > j. \qquad (114)$$

Let

$$\delta = [\alpha_1 / \|\alpha_1\|, \alpha_2 / \|\alpha_2\|, \dots, \alpha_r / \|\alpha_r\|], \qquad (115)$$

and

$$V = \begin{bmatrix} \|\alpha_1\| & \|\alpha_1\| u_{12} & \|\alpha_1\| u_{13} & \dots & \|\alpha_1\| u_{1r} \\ 0 & \|\alpha_2\| & \|\alpha_2\| u_{23} & \dots & \|\alpha_2\| u_{2r} \\ 0 & 0 & \|\alpha_3\| & \dots & \|\alpha_3\| u_{3r} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \|\alpha_r\| \end{bmatrix}. \qquad (116)$$

Then,

$$A_r = \delta V. \qquad (117)$$

Since $\delta$ is an orthonormal matrix, we have

$$\delta^T = \delta^{-1}.$$

Therefore,

$$\delta^T \delta = I \quad \text{and} \quad \delta \delta^T = \text{symmetric.}$$

In addition, since $V$ is an upper-triangular and nonsingular matrix, its ordinary inverse exists. Hence,

$$V^{-1} V = I.$$

It can be seen then that

$$A_r^+ = V^{-1} \delta^T, \qquad (118)$$
because the four properties of the Moore-Penrose generalized inverse are fulfilled. That is,

$$A_r^+ A_r = V^{-1} \delta^T \delta V = V^{-1} V = I, \qquad (119)$$

$$A_r A_r^+ = \delta V V^{-1} \delta^T = \delta \delta^T = \text{symmetric}, \qquad (120)$$

$$A_r^+ A_r A_r^+ = I A_r^+ = A_r^+, \qquad (121)$$

$$A_r A_r^+ A_r = A_r I = A_r. \qquad (122)$$

It follows that

$$x^{(r)} = A_r^+ d = V^{-1} \delta^T d, \qquad (123)$$

and furthermore,

$$g_r(d_r) = x^{(r)T} x^{(r)} = d^T \delta (V^T V)^{-1} \delta^T d. \qquad (124)$$

Comparing Eq. (124) against Eq. (58), we see that

$$R_r = \delta(V^T V)^{-1} \delta^T \qquad (125)$$

is an alternative formula for calculating $R_r$. That is to say, after carrying out the $\alpha Q$-algorithm, the value of $R_r$ can be directly obtained from the vectors $\alpha_1, \alpha_2, \dots, \alpha_r$ and $\alpha_1, \alpha_2, \dots, \alpha_r$ without considering the matrix inversion problem. In this sense, Eq. (125) has advantages over Eq. (57).

Furthermore, in view of the fact that the $\alpha Q \beta R$-algorithm solves the problem of finding the shortest length vector $x$ such that $\|Ax-b\|^2$ is a minimum, we can use it to determine the $j$th column of $A^+$, where $A$ is the matrix with its first $r$ columns being independent. This is because, when we put

$$b = [0, \dots, 0, 1, 0, \dots, 0]^T,$$

with the one in the $j$th position, the shortest length solution to

$$\|Ax-b\|^2 = \min,$$

if we call it $x^{(j)opt}$, is the $j$th column of $A^+$, since

$$x^{(j)opt} = A^+[0, \dots, 0, 1, 0, \dots, 0]^T.$$

This way, we can obtain all the columns of $A^+$ by letting $j=1, 2, \dots, m$. For a more general matrix $B$ with rank $r$, we may first use the $\alpha Q$ algorithm to transform it into the form of $A$, so that

$$B = AP,$$
where $\mathbf{P}$ is a permutation matrix. Since

$$ \mathbf{B}^{+} = (\mathbf{A}\mathbf{P})^{+} = \mathbf{P}^{+}\mathbf{A}^{+} = \mathbf{P}^{T}\mathbf{A}^{+}, $$

we see that the generalized inverse of $\mathbf{B}$ can be readily obtained by multiplying the matrix $\mathbf{A}^{+}$ by $\mathbf{P}^{T}$ on the left. Thus, the $\alpha\mathbf{Q}\beta\mathbf{R}$-algorithm is also an algorithm for calculating the generalized inverse of any matrix.

## 7. Concluding Remarks

An algorithm based on dynamic programming has now been given for finding shortest length solutions of least-squares problems in which the matrix may be of less than full rank. It provides a sequential solution to the problem. Numerical experiments need to be carried out to test the algorithm and to compare it against alternative procedures.

## References

1. KALABA, R., and RASAKHOO, N., *Algorithms for Generalized Inverses*, Journal of Optimization Theory and Applications, Vol. 48, No. 3, pp. 427-435, 1986.

2. LAWSON, C., and HANSON, R., *Solving Least-Squares Problems*, Prentice-Hall, Englewood Cliffs, New Jersey, 1974.

3. BELLMAN, R., *Dynamic Programming*, Princeton University Press, Princeton, New Jersey, 1957.

4. BELLMAN, R., KALABA, R., and LOCKETT, J., *Numerical Inversion of Laplace Transforms*, Elsevier, New York, New York, 1964.