Numerical strategies for recursive least squares solutions to the matrix equation AX = B

Abstract

The recursive solution to the Procrustes problem -with or without constraints- is thoroughly investigated. Given known matrices $\mathit{A}$ and $\mathit{B}$, the proposed solution minimizes the square of the Frobenius norm of the difference $\mathit{A}\mathit{X}-\mathit{B}$ when rows or columns are added to $\mathit{A}$ and $\mathit{B}$. The proposed method is based on efficient strategies which reduce the computational cost by utilizing previous computations when new data are acquired. This is particularly useful in the iterative solution of an unbalanced orthogonal Procrustes problem. The results show that the computational efficiency of the proposed recursive algorithms is more significant when the dimensions of the matrices are large. This demonstrates the usefulness of the proposed algorithms in the presence of high-dimensional data sets. The practicality of the new method is demonstrated through an application in machine learning, namely feature extraction for image processing.

Keywords:

• Orthogonal matrix
• large-scale matrix
• sequential updating
• high-dimensional data

MSC2020 AMS SUBJECT CLASSIFICATIONS:

• 62L12
• 15B10
• 15A23
• 68U10
• 15A24

1. Introduction

In practice, one may be interested in finding the matrix $\mathit{X}$ such that $\mathit{A}\mathit{X}=\mathit{B}$ where matrices $\mathit{A}$ and $\mathit{B}$ come from experiments. However, $\mathit{A}$ and $\mathit{B}$ often do not satisfy the solvability conditions and hence, the least squares solution of the difference $\mathit{A}\mathit{X}-\mathit{B}$ is required [15]. Specifically, the problem of approximating one given matrix $\mathit{A}$ with another given matrix $\mathit{B}$ by a transformation matrix $\mathit{X}$ so that the square of the difference $\mathit{A}\mathit{X}-\mathit{B}$ is minimized is known as the Procrustes problem [13,15,33]. Often, depending on the application, it is assumed that $\mathit{X}$ belongs to a specific class of matrices, and thus setting in this way a set of constraints to the optimization problem. The most frequent classes of matrices for $\mathit{X}$ is orthogonality and symmetry, and variants thereof, see, for example [5,11,15,19,21,29,31,33]. In many cases, orthogonal factorizations like the singular value decomposition (SVD), the eigenvalue decomposition (EVD) and the CS decomposition have been used to solve the Procrustes problem and variants thereof [8,15,19,22,31,33, pp. 327–328].

The application of the Procrustes problem in factor analysis has a long history [10,13,15,28]. It also appears in numerical analysis problems for the solution of partial differential equations, in multidimensional scaling, in growth curve modelling, in scientific computing, in computer vision, in image processing, in system and control theory, in the analysis of space structures and in aerospace engineering for spacecraft attitude determination [12,23,24,26,27,32,38,41,43]. Moreover, in the case where $\mathit{A}$ is the identical matrix the problem becomes a matrix nearness problem with many applications in statistical and financial modelling and in theoretical computer science, see, for example [20,36,39].

The recursive solution to a least squares problem is needed when the experiment is conducted repeatedly and as a result the given matrices are updated with new arriving data regularly. Also, in high dimensional settings, the matrices $\mathit{A}$ and $\mathit{B}$ are very large and it may not be possible to treat all data at once or the computational cost of processing them may be significantly expensive. In this case, a sequential procedure which splits $\mathit{A}$ and $\mathit{B}$ into sub-matrices of smaller dimensions and then proceeds by gradually incorporating the sub-matrices into the least squares solution of $\mathit{A}\mathit{X}-\mathit{B}$ is essential. Recursive least squares is often needed in many problems of different areas like engineering, statistics, econometrics and finance [7,9,17,18,42]. A recursive algorithm reduces the computational cost and also the storage requirements for large matrices.

Herein, the recursive least squares solution to the matrix equation $\mathit{A}\mathit{X}=\mathit{B}$ when $\mathit{A}$ and $\mathit{B}$ are known matrices is investigated in depth. Namely, the recursive solution to the Procrustes problem is examined. The use of the QR decomposition is examined when there are no constraints on $\mathit{X}$. Also, the problems of minimizing the difference $\mathit{A}\mathit{X}-\mathit{B}$ when $\mathit{X}$ is orthogonal and when $\mathit{X}$ is symmetric are also considered and their recursive least squares solution using the eigenvalue and singular value decompositions is explored in depth. When constraints are imposed on $\mathit{X}$, the method of Lagrange multipliers is used to solve the optimization problem. The proposed solution, in each case, is the matrix which minimizes the square of the Frobenious norm of $\mathit{A}\mathit{X}-\mathit{B}$. It is an exact solution in the sense that $\mathit{X}$ is explicitly determined and does not comprise any arbitrary elements. The recursive numerical solution proposed does not require the matrices be full rank.

Throughout this paper, $\Vert \cdot {\Vert }_F$ denotes the Frobenius norm. Also, for known matrices $\mathit{S}$ and $\mathit{P}$, when computing partial derivatives [25, p. 201] the following properties are used: (1) $\frac{\mathrm{\partial }\left(\mathit{S}\mathit{X}\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}={\mathit{S}}^T{\mathit{P}}^T,\frac{\mathrm{\partial }\left(\mathit{S}{\mathit{X}}^T\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}=\mathit{P}\mathit{S}.$(1) The paper is organized as follows. Section 2 introduces the general Procrustes problem where no assumption is made for the solution matrix $\mathit{X}$. The problem is solved using the QR decomposition and then the recursive solution is presented. Section 3 considers the orthogonal Procrustes problem where the solution matrix $\mathit{X}$ is orthogonal and Section 4 derives the solution to the symmetric Procrustes problem when $\mathit{X}$ is assumed to be symmetric. Section 5 presents computational results and finally, in Section 6 we conclude and discuss future work.

2. Numerical solution to the general procrustes problem

Consider the problem of finding a matrix $\mathit{X}\in {\mathbb{R}}^{n\times n}$ so that the known matrix $\mathit{B}\in {\mathbb{R}}^{m\times n}$ is approximated by matrix $\mathit{A}\mathit{X}$ where $\mathit{A}\in {\mathbb{R}}^{m\times n}$ is also known. That is, a solution to the matrix equation (2) $\mathit{A}\mathit{X}=\mathit{B}$(2) is required. The least squares approximation problem to be solved is given by (3) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2,$(3) where $\begin{array}{rl}f\left(\mathit{X}\right)={\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2& =\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\right(\mathit{A}\mathit{X}-\mathit{B}{)}^T(\mathit{A}\mathit{X}-\mathit{B}))\\ & =\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}({\mathit{X}}^T{\mathit{A}}^T\mathit{A}\mathit{X}-{\mathit{X}}^T{\mathit{A}}^T\mathit{B}-{\mathit{B}}^T\mathit{A}\mathit{X}+{\mathit{B}}^T\mathit{B}).\end{array}$On using (1(1) $\frac{\mathrm{\partial }\left(\mathit{S}\mathit{X}\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}={\mathit{S}}^T{\mathit{P}}^T,\frac{\mathrm{\partial }\left(\mathit{S}{\mathit{X}}^T\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}=\mathit{P}\mathit{S}.$(1) ) partial differentiation yields (4) $\frac{\mathrm{\partial }f\left(\mathit{X}\right)}{\mathrm{\partial }\mathit{X}}=2{\mathit{A}}^T\mathit{A}\mathit{X}-2{\mathit{A}}^T\mathit{B}=\mathbf{0}$(4) whence (5) ${\mathit{A}}^T\mathit{A}\mathit{X}={\mathit{A}}^T\mathit{B}.$(5) When ${\mathit{A}}^T\mathit{A}$ is non-singular the solution to (5(5) ${\mathit{A}}^T\mathit{A}\mathit{X}={\mathit{A}}^T\mathit{B}.$(5) ) is given by $\mathit{X}=({\mathit{A}}^T\mathit{A}{)}^{-1}{\mathit{A}}^T\mathit{B}$ [14]. However, ${\mathit{A}}^T\mathit{A}$ may be ill-conditioned and in this case inverting the matrix may give inaccurate results. In the case where ${\mathit{A}}^T\mathit{A}$ is singular the latter will fail to give a solution for $\mathit{X}$. A numerically stable method to obtain $\mathit{X}$ is to use the QR decomposition (QRD) of $\mathit{A}$, namely (6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) Then $\mathit{A}={\mathit{Q}}_{A1}{\mathit{R}}_A$ and hence ${\mathit{A}}^T\mathit{A}={\mathit{R}}_A^T{\mathit{R}}_A$ and ${\mathit{A}}^T\mathit{B}={\mathit{R}}_A^T{\mathit{R}}_{B1}$. Therefore, $\mathit{X}={\mathit{R}}_A^{-1}{\mathit{R}}_{B1}$. When $\mathit{A}$ is rank deficient, the procedure to obtain $\mathit{X}$ is similar to the above, namely (6(6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) ). In this case, a complete QRD can be computed to triangularize $\mathit{A}$ [1].

2.1. Recursive solution

We next consider the recursive solution of the Procrustes problem when the matrices $\mathit{A}$ and $\mathit{B}$ are updated. Without loss of generality, suppose that $\mathit{A}$ and $\mathit{B}$ are augmented with the addition of a single row; namely, (7) $\stackrel{\mathbf{~}}{\mathit{A}}=\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\stackrel{\mathbf{~}}{\mathit{B}}=\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right),$(7) where $\mathit{A},\mathit{B}$ are as in (2(2) $\mathit{A}\mathit{X}=\mathit{B}$(2) ) and $\mathit{a},\mathit{b}\in {\mathbb{R}}^{1\times n}$ represent new data points. Then, the updated Procrustes problem, based on $\stackrel{\mathbf{~}}{\mathit{A}}$ and $\stackrel{\mathbf{~}}{\mathit{B}}$, requires the solution of the least squares problem (8) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2$(8) when (3(3) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2,$(3) ) has already been solved (see (4(4) $\frac{\mathrm{\partial }f\left(\mathit{X}\right)}{\mathrm{\partial }\mathit{X}}=2{\mathit{A}}^T\mathit{A}\mathit{X}-2{\mathit{A}}^T\mathit{B}=\mathbf{0}$(4) )–(6(6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) )). The efficient solution of (8(8) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2$(8) ) requires that previous computations from the solution of (3(3) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2,$(3) ) be utilized. Namely, $\begin{array}{rl}\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2& =\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \left(\begin{array}{cc}{\mathit{Q}}_{A1}^T& \mathbf{0}\\ \mathbf{0}& {\mathbf{1}}_n\end{array}\right)[\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\tilde{\mathit{X}}-\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right)]\Vert }_F^2\\ & =\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \left(\begin{array}{c}{\mathit{R}}_A\\ \mathit{a}\end{array}\right)\tilde{\mathit{X}}-\left(\begin{array}{c}{\mathit{R}}_{B1}\\ \mathit{b}\end{array}\right)\Vert }_F^2,\end{array}$where ${\mathit{Q}}_{A1}$, ${\mathit{R}}_A$ and ${\mathit{R}}_{B1}$ are as in (6(6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) ) and ${\mathbf{1}}_n$ is the n-dimensional row vector of ones. Consider now the updating QRD (9) ${\mathit{Q}}_{A_u}^T\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathit{a}& \mathit{b}\end{array}\right)=\left(\begin{array}{cc}{\tilde{\mathit{R}}}_A& {\tilde{\mathit{R}}}_{B_1}\\ \mathbf{0}& {\tilde{\mathit{R}}}_{B_2}\end{array}\right).$(9) We then have that $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2=\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}({\Vert {\stackrel{\mathbf{~}}{\mathit{R}}}_A\stackrel{\mathbf{~}}{\mathit{X}}-{\stackrel{\mathbf{~}}{\mathit{R}}}_{B_1}\Vert }_F^2+{\Vert {\stackrel{\mathbf{~}}{\mathit{R}}}_{B_2}\Vert }_F^2)$and therefore it follows that $\stackrel{\mathbf{~}}{\mathit{X}}={\stackrel{\mathbf{~}}{\mathit{R}}}_A^{-1}{\stackrel{\mathbf{~}}{\mathit{R}}}_{B_1}$.

In many cases, it is possible that the matrices $\mathit{A}$ and $\mathit{B}$ are updated with a single column, namely $\stackrel{\mathbf{ˇ}}{\mathit{A}}=\left(\begin{array}{cc}\mathit{A}& {\mathit{A}}_{n+1}\end{array}\right),\stackrel{\mathbf{ˇ}}{\mathit{B}}=\left(\begin{array}{cc}\mathit{B}& {\mathit{B}}_{n+1}\end{array}\right)$where ${\mathit{A}}_{n+1},{\mathit{B}}_{n+1}\in {\mathbb{R}}^{m\times 1}$ denote new variables which become available after (3(3) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2,$(3) ) has been solved. As a result, the solution $\stackrel{\mathbf{ˇ}}{\mathit{X}}\in {\mathbb{R}}^{(n+1)\times (n+1)}$ to the updated Procrustes problem (10) $\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2.$(10) needs to be computed. By utilizing efficiently previous computations from the solution of (3(3) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2,$(3) ), (10(10) $\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2.$(10) ) is written as $\begin{array}{rl}\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2& =\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert {\mathit{Q}}_A^T(\left(\begin{array}{cc}\mathit{A}& {\mathit{A}}_{n+1}\end{array}\right)\stackrel{\mathbf{ˇ}}{\mathit{X}}-\left(\begin{array}{cc}\mathit{B}& {\mathit{B}}_{n+1}\end{array}\right))\Vert }_F^2\\ & =\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \left(\begin{array}{cc}{\mathit{R}}_A& {\stackrel{\mathbf{~}}{\mathit{A}}}_{n+1}\\ 0& {\stackrel{\mathbf{ˆ}}{\mathit{A}}}_{n+1}\end{array}\right)\stackrel{\mathbf{ˇ}}{\mathit{X}}-\left(\begin{array}{cc}{\mathit{R}}_{B_1}& {\stackrel{\mathbf{~}}{\mathit{B}}}_{n+1}\\ {\mathit{R}}_{B_2}& {\stackrel{\mathbf{ˆ}}{\mathit{B}}}_{n+1}\end{array}\right)\Vert }_F^2,\end{array}$where ${\mathit{Q}}_A$ is from the QRD of $\mathit{A}$ in (6(6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) ). The column-updating QRD that needs to be computed is then given by $\left(\begin{array}{cc}{\mathit{I}}_n& \mathbf{0}\\ \mathbf{0}& {\stackrel{\mathbf{ˇ}}{\mathit{q}}}^T\end{array}\right)\left(\begin{array}{cc}{\mathit{R}}_A& {\stackrel{\mathbf{~}}{\mathit{A}}}_{n+1}\\ 0& {\stackrel{\mathbf{ˆ}}{\mathit{A}}}_{n+1}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\stackrel{\mathbf{~}}{\mathit{A}}}_{n+1}\\ 0& \check{a}\\ 0& 0\end{array}\right),$where $\stackrel{\mathbf{ˇ}}{\mathit{q}}$ is an orthogonal transformation that eliminates all but the first element of ${\stackrel{\mathbf{ˆ}}{\mathit{A}}}_{n+1}$ and $\check{a}$ is a scalar. Hence, the updated Procrustes problem (10(10) $\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2.$(10) ) becomes (11) $\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2=\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}({\Vert {\stackrel{\mathbf{ˇ}}{\mathit{R}}}_A\stackrel{\mathbf{ˇ}}{\mathit{X}}-{\stackrel{\mathbf{ˇ}}{\mathit{R}}}_B\Vert }_F^2+{\Vert {\stackrel{\mathbf{ˇ}}{\mathit{R}}}_{B_2}\Vert }_F^2),$(11) where ${\stackrel{\mathbf{ˇ}}{\mathit{R}}}_A=\left(\begin{array}{cc}{\mathit{R}}_A& {\stackrel{\mathbf{~}}{\mathit{A}}}_{n+1}\\ 0& \check{a}\end{array}\right),{\stackrel{\mathbf{ˇ}}{\mathit{R}}}_B=\left(\begin{array}{cc}{\mathit{R}}_{B_1}& {\stackrel{\mathbf{~}}{\mathit{B}}}_{n+1}\\ {\mathit{R}}_{B_2}^{\left(1\right)}& {\stackrel{\mathbf{ˆ}}{\mathit{B}}}_{n+1}^{\left(1\right)}\end{array}\right)\mathrm{a}\mathrm{n}\mathrm{d}{\stackrel{\mathbf{ˇ}}{\mathit{R}}}_{B_2}=\left(\begin{array}{cc}{\mathit{R}}_{B_2}^{\left(2\right)}& {\stackrel{\mathbf{ˆ}}{\mathit{B}}}_{n+1}^{\left(2\right)}\end{array}\right).$The solution to problem (11(11) $\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{ˇ}}{\mathit{A}}\stackrel{\mathbf{ˇ}}{\mathit{X}}-\stackrel{\mathbf{ˇ}}{\mathit{B}}\Vert }_F^2=\underset{\stackrel{\mathbf{ˇ}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}({\Vert {\stackrel{\mathbf{ˇ}}{\mathit{R}}}_A\stackrel{\mathbf{ˇ}}{\mathit{X}}-{\stackrel{\mathbf{ˇ}}{\mathit{R}}}_B\Vert }_F^2+{\Vert {\stackrel{\mathbf{ˇ}}{\mathit{R}}}_{B_2}\Vert }_F^2),$(11) ) is given by $\stackrel{\mathbf{ˇ}}{\mathit{X}}={\stackrel{\mathbf{ˇ}}{\mathit{R}}}_A^{-1}{\stackrel{\mathbf{ˇ}}{\mathit{R}}}_B$.

3. The orthogonal procrustes problem

The orthogonal Procrustes problem (OPP) is that of minimizing the sum of the squared error of the difference matrix $\mathit{A}\mathit{X}-\mathit{B}$ when the unknown matrix $\mathit{X}$ is orthogonal. The constraint is imposed by using the method of Lagrange multipliers; the matrices $\mathit{A}$ and $\mathit{B}$ need not be full rank [33]. The constrained optimization problem is then given by (12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) where $\mathit{A},\mathit{B}\in {\mathbb{R}}^{m\times n}$ and $\mathit{X}\in {\mathbb{R}}^{n\times n}$ is orthogonal. Herein, we are most interested in cases where m<n. To find the solution to (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ), we consider the Lagrangian function $\mathcal{L}\left(\mathit{X}\right)=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\right(\mathit{A}\mathit{X}-\mathit{B}{)}^T(\mathit{A}\mathit{X}-\mathit{B}))+\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\mathbf{\Lambda }\right(\mathit{X}{\mathit{X}}^T-\mathit{I}\left)\right),$which is equivalently written as (13) $\mathcal{L}\left(\mathit{X}\right)=\sum _{j=1}^n\sum _{i=1}^m{(\sum _{q=1}^n{a}_{iq}{x}_{qj}-b_{ij})}^2+\sum _{q,r=1}^n{\lambda }_{qr}(\sum _{k=1}^n{x}_{qs}{x}_{rs}-{\delta }_{qr}),$(13) where $\mathbf{\Lambda }=[{\lambda }_{qr}{]}_{q,r=1}^n$ is the symmetric matrix of Lagrange multipliers [15]. Partial differentiation of (13(13) $\mathcal{L}\left(\mathit{X}\right)=\sum _{j=1}^n\sum _{i=1}^m{(\sum _{q=1}^n{a}_{iq}{x}_{qj}-b_{ij})}^2+\sum _{q,r=1}^n{\lambda }_{qr}(\sum _{k=1}^n{x}_{qs}{x}_{rs}-{\delta }_{qr}),$(13) ) yields (14) $\frac{\mathrm{\partial }\mathcal{L}\left(\mathit{X}\right)}{\mathrm{\partial }{x}_{pj}}=2\sum _{i=1}^m(\sum _{q=1}^n{a}_{iq}{x}_{qj}-b_{ij})a_{ip}+2\sum _{r=1}^n{\lambda }_{qr}{x}_{rj}.$(14) On setting (14(14) $\frac{\mathrm{\partial }\mathcal{L}\left(\mathit{X}\right)}{\mathrm{\partial }{x}_{pj}}=2\sum _{i=1}^m(\sum _{q=1}^n{a}_{iq}{x}_{qj}-b_{ij})a_{ip}+2\sum _{r=1}^n{\lambda }_{qr}{x}_{rj}.$(14) ) equal to zero and using matrix notation, (14(14) $\frac{\mathrm{\partial }\mathcal{L}\left(\mathit{X}\right)}{\mathrm{\partial }{x}_{pj}}=2\sum _{i=1}^m(\sum _{q=1}^n{a}_{iq}{x}_{qj}-b_{ij})a_{ip}+2\sum _{r=1}^n{\lambda }_{qr}{x}_{rj}.$(14) ) can be written as ${\mathit{A}}_q^T{\mathit{A}}_q{\mathit{X}}_j+{\mathbf{\Lambda }}_q{\mathit{X}}_j={\mathit{A}}_q^T{\mathit{B}}_j,q,j=1,\dots ,n,$or equivalently, as (15) $({\mathit{A}}^T\mathit{A}+\mathbf{\Lambda })\mathit{X}={\mathit{A}}^T\mathit{B},$(15) where ${\mathit{A}}_q$ is the qth column of $\mathit{A}$, ${\mathit{X}}_j$, ${\mathit{B}}_j$ are the jth columns of $\mathit{X}$ and $\mathit{B}$, respectively, and ${\mathbf{\Lambda }}_q$ is the qth row of $\mathbf{\Lambda }$. From (15(15) $({\mathit{A}}^T\mathit{A}+\mathbf{\Lambda })\mathit{X}={\mathit{A}}^T\mathit{B},$(15) ), $({\mathit{A}}^T\mathit{A}+\mathbf{\Lambda })\mathit{X}{\mathit{X}}^T({\mathit{A}}^T\mathit{A}+\mathbf{\Lambda }{)}^T={\mathit{A}}^T\mathit{B}({\mathit{A}}^T\mathit{B}{)}^T$, therefore $\mathbf{\Lambda }=({\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}{)}^{1/2}-{\mathit{A}}^T\mathit{A}$ and, as observed earlier, $\mathbf{\Lambda }$ is symmetric. Now on post-multiplying (15(15) $({\mathit{A}}^T\mathit{A}+\mathbf{\Lambda })\mathit{X}={\mathit{A}}^T\mathit{B},$(15) ) by ${\mathit{X}}^T$ gives ${\mathit{A}}^T\mathit{A}+\mathbf{\Lambda }={\mathit{A}}^T\mathit{B}{\mathit{X}}^T$ which implies that $\mathbf{\Lambda }={\mathit{A}}^T{\mathit{B}\mathit{X}}^T-{\mathit{A}}^T\mathit{A}$, whence ${\mathit{A}}^T{\mathit{B}\mathit{X}}^T=\mathit{X}{\mathit{B}}^T\mathit{A}$. Therefore, (16) ${\mathit{A}}^T\mathit{B}=\mathit{X}{\mathit{B}}^T\mathit{A}\mathit{X}.$(16) Furthermore, let $\mathit{H}={\mathit{A}}^T\mathit{B}$ and consider the following two matrices $\mathit{F}=\mathit{H}{\mathit{H}}^T={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}$and $\mathit{G}={\mathit{H}}^T\mathit{H}={\mathit{B}}^T\mathit{A}{\mathit{A}}^T\mathit{B}.$Matrices $\mathit{F},\mathit{G}\in {\mathbb{R}}^{n\times n}$ are symmetric and thus they are diagonalizable and their eigenvalue decomposition (EVD) exists, that is, (17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) and (17b) $\mathit{G}=\mathit{V}\mathit{D}{\mathit{V}}^T,$(17b) where $\mathit{U},\mathit{V}\in {\mathbb{R}}^{n\times n}$ are orthogonal. Additionally, since they are of the form $\mathit{H}{\mathit{H}}^T$ and ${\mathit{H}}^T\mathit{H}$ they have the same eigenvalues. Now on using (16(16) ${\mathit{A}}^T\mathit{B}=\mathit{X}{\mathit{B}}^T\mathit{A}\mathit{X}.$(16) ) it follows that $\begin{array}{rl}\mathit{F}& ={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}\\ & =\left(\mathit{X}{\mathit{B}}^T\mathit{A}\mathit{X}\right){\left(\mathit{X}{\mathit{B}}^T\mathit{A}\mathit{X}\right)}^T\\ & =\mathit{X}{\mathit{B}}^T\mathit{A}{\mathit{A}}^T\mathit{B}{\mathit{X}}^T\\ & =\mathit{X}\mathit{G}{\mathit{X}}^T,\end{array}$and from (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ) and (17b(17b) $\mathit{G}=\mathit{V}\mathit{D}{\mathit{V}}^T,$(17b) ) we now have that $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T=\mathit{X}\mathit{V}\mathit{D}{\mathit{V}}^T{\mathit{X}}^T,$which implies that $\mathit{U}=\mathit{X}\mathit{V}$, where $\mathit{U}$ is as in (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ). Therefore, the solution to the least squares problem (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) is given by $\mathit{X}=\mathit{U}{\mathit{V}}^T.$Finally, a sufficient condition which guaranties the uniqueness of the solution $\mathit{X}$ and that the argument in (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) is minimized requires that all the diagonal elements of the matrix ${\mathit{D}}^{1/2}$ in (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ) and (17b(17b) $\mathit{G}=\mathit{V}\mathit{D}{\mathit{V}}^T,$(17b) ) are non-negative [15,33]. In essence, this condition determines the orientation of the orthogonal matrices $\mathit{U}$ and $\mathit{V}$ and is specified by the Eckart–Young decomposition (see [10]) of ${\mathit{A}}^T\mathit{B}$, namely ${\mathit{A}}^T\mathit{B}=\mathit{U}{\mathit{D}}^{1/2}{\mathit{V}}^T.$Notice that in the case of symmetric orthogonality the procedure is the same but at the last step the symmetry of $\mathit{X}$ needs to be taken into account. Namely, the sum of $\mathit{U}{\mathit{V}}^T$ is not explicitly computed since only the upper or lower part of $\mathit{X}$ needs to be determined. This results in a dimensional reduction of the solution matrix from $n^2$ to $n(n+1)/2$ [40].

3.1. Recursive solution to the orthogonal procrustes problem

Suppose that an updated orthogonal Procrustes problem needs to be solved when new data become available. Without loss of generality it is assumed that the original OPP (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) needs to be re-solved when appending a new single row of data in matrices $\mathit{A}$ and $\mathit{B}$. That is, let the row updated matrices $\stackrel{\mathbf{~}}{\mathit{A}}$ and $\stackrel{\mathbf{~}}{\mathit{B}}$ be as in (7(7) $\stackrel{\mathbf{~}}{\mathit{A}}=\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\stackrel{\mathbf{~}}{\mathit{B}}=\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right),$(7) ), where $\mathit{A},\mathit{B}$ are as in (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) with $\mathit{a},\mathit{b}\in {\mathbb{R}}^{1\times n}$. The updated orthogonal Procrustes problem based on $\stackrel{\mathbf{~}}{\mathit{A}}$ and $\stackrel{\mathbf{~}}{\mathit{B}}$ requires that the solution of the following least squares problem be derived: (18) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}\stackrel{\mathbf{~}}{\mathit{X}}}^T={\stackrel{\mathbf{~}}{\mathit{X}}}^T\stackrel{\mathbf{~}}{\mathit{X}}=\mathit{I},$(18) where $\stackrel{\mathbf{~}}{\mathit{X}}\in {\mathbb{R}}^{n\times n}$. The solution to the least squares problem (18(18) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}\stackrel{\mathbf{~}}{\mathit{X}}}^T={\stackrel{\mathbf{~}}{\mathit{X}}}^T\stackrel{\mathbf{~}}{\mathit{X}}=\mathit{I},$(18) ) is obtained by using the method of Lagrange multipliers, namely (19) $\mathcal{L}\left(\stackrel{\mathbf{~}}{\mathit{X}}\right)=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}(\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}{)}^T(\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\left)\right)+\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\stackrel{\mathbf{~}}{\mathbf{\Lambda }}\right(\stackrel{\mathbf{~}}{\mathit{X}}{\stackrel{\mathbf{~}}{\mathit{X}}}^T-\mathit{I}\left)\right)$(19) which yields the first order condition $\frac{\mathrm{\partial }\mathcal{L}\left(\stackrel{\mathbf{~}}{\mathit{X}}\right)}{\mathrm{\partial }\stackrel{\mathbf{~}}{\mathit{X}}}=2{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-2{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}+({\stackrel{\mathbf{~}}{\mathbf{\Lambda }}}^T+\stackrel{\mathbf{~}}{\mathbf{\Lambda }})\stackrel{\mathbf{~}}{\mathit{X}}$on using (1(1) $\frac{\mathrm{\partial }\left(\mathit{S}\mathit{X}\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}={\mathit{S}}^T{\mathit{P}}^T,\frac{\mathrm{\partial }\left(\mathit{S}{\mathit{X}}^T\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}=\mathit{P}\mathit{S}.$(1) ). In a way similar to (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ), (19(19) $\mathcal{L}\left(\stackrel{\mathbf{~}}{\mathit{X}}\right)=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}(\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}{)}^T(\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\left)\right)+\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\stackrel{\mathbf{~}}{\mathbf{\Lambda }}\right(\stackrel{\mathbf{~}}{\mathit{X}}{\stackrel{\mathbf{~}}{\mathit{X}}}^T-\mathit{I}\left)\right)$(19) ) has the solution (20) $\stackrel{\mathbf{~}}{\mathit{X}}=\stackrel{\mathbf{~}}{\mathit{U}}{\stackrel{\mathbf{~}}{\mathit{V}}}^T,$(20) where $\stackrel{\mathbf{~}}{\mathit{U}}$ and $\stackrel{\mathbf{~}}{\mathit{V}}$ are the orthogonal matrices of the EVD of ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$ and ${\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}$, respectively.

The recursive solution of the OPP presumes that previous computations in (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ) from the solution of the original Procrustes problem (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) are efficiently utilized. Consider the recursive computation of the EVD of ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$, namely (21) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& =({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b})({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b}{)}^T\\ & ={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\end{array}$(21) on using (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ). Therefore, the recursive solution of an orthogonal Procrustes problem becomes a modified symmetric matrix eigenvalue problem [3,4,16]. In particular, (21(21) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& =({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b})({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b}{)}^T\\ & ={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\end{array}$(21) ) implies three rank-1 modifications of the EVD ${\mathit{U}\mathit{D}\mathit{U}}^T$ of the matrix ${\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}$ in (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ). That is, the EVD of ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$ is obtained recursively in three main steps. First, (21(21) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& =({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b})({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b}{)}^T\\ & ={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\end{array}$(21) ) is written as (22) ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\left(\begin{array}{cc}0& 1\\ 1& \text{β}\end{array}\right)\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right),$(22) where $\tilde{\mathit{b}}=\mathit{b}{\mathit{B}}^{\mathit{T}}\mathit{A}$ and $\text{β}={\mathit{b}\mathit{b}}^T$ is a scalar since $\mathit{b}$ and $\mathit{a}$ are row vectors. Second, the following EVD is derived (23) $\mathbf{\Delta }=\left(\begin{array}{cc}0& 1\\ 1& \text{β}\end{array}\right)=\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T$(23) and, since $\mathbf{\Delta }$ is a symmetric matrix, it is diagonalizable. Third, on using  (22(22) ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\left(\begin{array}{cc}0& 1\\ 1& \text{β}\end{array}\right)\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right),$(22) ), (23(23) $\mathbf{\Delta }=\left(\begin{array}{cc}0& 1\\ 1& \text{β}\end{array}\right)=\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T$(23) ) becomes (24) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)}^T\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right)\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\text{β}}_1{\tilde{\mathit{b}}}_1^T{\tilde{\mathit{b}}}_1+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2\\ & =\mathit{U}(\mathit{D}+{\text{β}}_1{\widetilde{\stackrel{~}{\mathit{b}}}}_1^T{\widetilde{\stackrel{~}{\mathit{b}}}}_1){\mathit{U}}^T+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2,\end{array}$(24) where ${\widetilde{\stackrel{~}{\mathit{b}}}}_1={\tilde{\mathit{b}}}_1\mathit{U}$ with ${\text{β}}_1$, ${\text{β}}_2$ being the eigenvalues of $\mathbf{\Delta }$. Consider now the sequential updating of the diagonal matrix $\mathit{D}$ in two steps. The first step computes the EVD of $\mathit{D}+{\text{β}}_1{\widetilde{\stackrel{~}{\mathit{b}}}}_1^T{\widetilde{\stackrel{~}{\mathit{b}}}}_1={\mathit{U}}_1{\mathit{D}}_1{\mathit{U}}_1^T$, where ${\mathit{U}}_1\in {\mathbb{R}}^{n\times n}$ is orthogonal and ${\mathit{D}}_1\in {\mathbb{R}}^{n\times n}$ is diagonal. The second step computes the EVD of ${\mathit{D}}_1+{\text{β}}_2{\widetilde{\stackrel{~}{\mathit{b}}}}_2^T{\widetilde{\stackrel{~}{\mathit{b}}}}_2={\mathit{U}}_2\tilde{\mathit{D}}{\mathit{U}}_2^T$, where ${\mathit{U}}_2\in {\mathbb{R}}^{n\times n}$ is orthogonal and $\tilde{\mathit{D}}\in {\mathbb{R}}^{n\times n}$ is the diagonal matrix with elements the eigenvalues of ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$. That is, ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}=\tilde{\mathit{U}}\tilde{\mathit{D}}{\tilde{\mathit{U}}}^T$, where $\tilde{\mathit{U}}=\mathit{U}{\mathit{U}}_1{\mathit{U}}_2$ and $\mathit{U}$ is as in (17b(17b) $\mathit{G}=\mathit{V}\mathit{D}{\mathit{V}}^T,$(17b) ). Repeating the procedure at steps (21(21) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& =({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b})({\mathit{A}}^T\mathit{B}+{\mathit{a}}^T\mathit{b}{)}^T\\ & ={\mathit{A}}^T\mathit{B}{\mathit{B}}^T\mathit{A}+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\mathit{A}}^T{\mathit{B}\mathit{b}}^T\mathit{a}+{\mathit{a}}^T{\mathit{b}\mathit{B}}^T\mathit{A}+{\mathit{a}}^T{\mathit{b}\mathit{b}}^T\mathit{a}\end{array}$(21) )–(24(24) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)}^T\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right)\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\text{β}}_1{\tilde{\mathit{b}}}_1^T{\tilde{\mathit{b}}}_1+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2\\ & =\mathit{U}(\mathit{D}+{\text{β}}_1{\widetilde{\stackrel{~}{\mathit{b}}}}_1^T{\widetilde{\stackrel{~}{\mathit{b}}}}_1){\mathit{U}}^T+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2,\end{array}$(24) ) for ${\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}$ will derive its EVD recursively and will therefore give ${\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}=\tilde{\mathit{V}}\tilde{\mathit{D}}{\tilde{\mathit{V}}}^T$, where $\tilde{\mathit{V}}=\mathit{V}{\mathit{V}}_1{\mathit{V}}_2$, $\mathit{V}$ is defined in (17a(17a) $\mathit{F}=\mathit{U}\mathit{D}{\mathit{U}}^T$(17a) ) and ${\mathit{V}}_1,{\mathit{V}}_2$ are computed in a way similar to that for (24(24) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)}^T\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right)\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\text{β}}_1{\tilde{\mathit{b}}}_1^T{\tilde{\mathit{b}}}_1+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2\\ & =\mathit{U}(\mathit{D}+{\text{β}}_1{\widetilde{\stackrel{~}{\mathit{b}}}}_1^T{\widetilde{\stackrel{~}{\mathit{b}}}}_1){\mathit{U}}^T+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2,\end{array}$(24) ). Therefore, the updated solution (20(20) $\stackrel{\mathbf{~}}{\mathit{X}}=\stackrel{\mathbf{~}}{\mathit{U}}{\stackrel{\mathbf{~}}{\mathit{V}}}^T,$(20) ) to the Procrustes problem (18(18) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}\stackrel{\mathbf{~}}{\mathit{X}}}^T={\stackrel{\mathbf{~}}{\mathit{X}}}^T\stackrel{\mathbf{~}}{\mathit{X}}=\mathit{I},$(18) ) has been derived.

When extra columns are added to $\mathit{A}$ and $\mathit{B}$, that is when (25) $\stackrel{\mathbf{˘}}{\mathit{A}}=\begin{array}{ccc}& n& k\\ & \mathit{A}& {\mathit{A}}_{n+1}\end{array}\mathrm{a}\mathrm{n}\mathrm{d}\stackrel{\mathbf{˘}}{\mathit{B}}=\begin{array}{ccc}& n& k\\ & \mathit{B}& {\mathit{B}}_{n+1}\end{array},$(25) an updated orthogonal Procrustes problem of larger dimensions needs to be solved. Namely, (26) $\underset{\stackrel{\mathbf{˘}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{˘}}{\mathit{A}}\stackrel{\mathbf{˘}}{\mathit{X}}-\stackrel{\mathbf{˘}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{˘}}{\mathit{X}}\stackrel{\mathbf{˘}}{\mathit{X}}}^T={\stackrel{\mathbf{˘}}{\mathit{X}}}^T\stackrel{\mathbf{˘}}{\mathit{X}}=\mathit{I},$(26) where $\stackrel{\mathbf{˘}}{\mathit{X}}\in {\mathbb{R}}^{(n+k)\times (n+k)}$ has been augmented by k columns and k rows. In a similar way as for (18(18) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}\stackrel{\mathbf{~}}{\mathit{X}}}^T={\stackrel{\mathbf{~}}{\mathit{X}}}^T\stackrel{\mathbf{~}}{\mathit{X}}=\mathit{I},$(18) ), the solution of (26(26) $\underset{\stackrel{\mathbf{˘}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{˘}}{\mathit{A}}\stackrel{\mathbf{˘}}{\mathit{X}}-\stackrel{\mathbf{˘}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{˘}}{\mathit{X}}\stackrel{\mathbf{˘}}{\mathit{X}}}^T={\stackrel{\mathbf{˘}}{\mathit{X}}}^T\stackrel{\mathbf{˘}}{\mathit{X}}=\mathit{I},$(26) ) is given by $\stackrel{\mathbf{˘}}{\mathit{X}}=\stackrel{\mathbf{˘}}{\mathit{U}}{\stackrel{\mathbf{˘}}{\mathit{V}}}^T$. It is obtained recursively, by updating the original orthogonal decompositions as in (24(24) $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)}^T\mathit{Q}\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right){\mathit{Q}}^T\left(\begin{array}{c}\tilde{\mathit{b}}\\ \mathit{a}\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)}^T\left(\begin{array}{cc}{\text{β}}_1& 0\\ 0& {\text{β}}_2\end{array}\right)\left(\begin{array}{c}{\tilde{\mathit{b}}}_1\\ {\tilde{\mathit{b}}}_2\end{array}\right)\\ & ={\mathit{U}\mathit{D}\mathit{U}}^T+{\text{β}}_1{\tilde{\mathit{b}}}_1^T{\tilde{\mathit{b}}}_1+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2\\ & =\mathit{U}(\mathit{D}+{\text{β}}_1{\widetilde{\stackrel{~}{\mathit{b}}}}_1^T{\widetilde{\stackrel{~}{\mathit{b}}}}_1){\mathit{U}}^T+{\text{β}}_2{\tilde{\mathit{b}}}_2^T{\tilde{\mathit{b}}}_2,\end{array}$(24) ). Notice that the addition of extra columns implies a rank-k updating of the EVD decomposition.

4. The symmetric procrustes problem

Consider the problem of minimizing the sum of squared error of the difference matrix ${\mathit{A}\mathit{X}}_S-\mathit{B}$ when the unknown matrix ${\mathit{X}}_S$ is symmetric. The constrained optimization problem is given by $\underset{{\mathit{X}}_S}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}{\mathit{X}}_S-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}}_S^T={\mathit{X}}_S,$where ${\mathit{X}}_S\in {\mathbb{R}}^{n\times n}$ is a symmetric matrix. As in the case of the OPP (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ), the matrices $\mathit{A}$ and $\mathit{B}$ are not necessarily full rank. Using the method of Lagrange multipliers, the problem becomes that of finding the matrix ${\mathit{X}}_S$ which minimizes $\mathcal{L}\left({\mathit{X}}_S\right)=\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\right(\mathit{A}{\mathit{X}}_{\mathit{S}}-\mathit{B}{)}^T(\mathit{A}{\mathit{X}}_{\mathit{S}}-\mathit{B}))+\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left(\mathbf{\Lambda }\right({\mathit{X}}_{\mathit{S}}-{\mathit{X}}_{\mathit{S}}^T\left)\right).$On using (1(1) $\frac{\mathrm{\partial }\left(\mathit{S}\mathit{X}\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}={\mathit{S}}^T{\mathit{P}}^T,\frac{\mathrm{\partial }\left(\mathit{S}{\mathit{X}}^T\mathit{P}\right)}{\mathrm{\partial }\mathit{X}}=\mathit{P}\mathit{S}.$(1) ) partial differentiation yields $\frac{\mathrm{\partial }L\left({\mathit{X}}_{\mathit{S}}\right)}{\mathrm{\partial }{\mathit{X}}_{\mathit{S}}}=2{\mathit{A}}^T\mathit{A}{\mathit{X}}_{\mathit{S}}-2{\mathit{A}}^T\mathit{B}+{\mathbf{\Lambda }}^T-\mathbf{\Lambda },$which is set to zero. Now let the matrix $\mathit{M}=\mathbf{\Lambda }-{\mathbf{\Lambda }}^T$, which is skew-symmetric, that is $\mathit{M}=-{\mathit{M}}^T$. It follows that $2{\mathit{A}}^T\mathit{A}{\mathit{X}}_{\mathit{S}}-2{\mathit{A}}^T\mathit{B}=-(2{\mathit{X}}_{\mathit{S}}{\mathit{A}}^T\mathit{A}-2{\mathit{B}}^T\mathit{A}),$which yields the Lyapunov equation (27) ${\mathit{A}}^T\mathit{A}{\mathit{X}}_{\mathit{S}}+{\mathit{X}}_{\mathit{S}}{\mathit{A}}^T\mathit{A}={\mathit{A}}^T\mathit{B}+{\mathit{B}}^T\mathit{A}.$(27) Since ${\mathit{A}}^T\mathit{A}$ is symmetric, it is therefore a diagonalizable matrix. That is, there is an orthogonal matrix $\mathit{P}$ and a diagonal matrix ${\mathit{D}}_A$ such that (28) ${\mathit{A}}^T\mathit{A}=\mathit{P}{\mathit{D}}_A{\mathit{P}}^T,$(28) where $\mathit{P}\in {\mathbb{R}}^{n\times n}$ has columns the eigenvectors of ${\mathit{A}}^{\mathit{T}}\mathit{A}$ and ${\mathit{D}}_A=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }_1,\dots ,{\mu }_n)$ has diagonal elements the eigenvalues of ${\mathit{A}}^{\mathit{T}}\mathit{A}$. Using  (27(27) ${\mathit{A}}^T\mathit{A}{\mathit{X}}_{\mathit{S}}+{\mathit{X}}_{\mathit{S}}{\mathit{A}}^T\mathit{A}={\mathit{A}}^T\mathit{B}+{\mathit{B}}^T\mathit{A}.$(27) ), (28(28) ${\mathit{A}}^T\mathit{A}=\mathit{P}{\mathit{D}}_A{\mathit{P}}^T,$(28) ) becomes (29) ${\mathit{D}}_A{\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}+{\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}{\mathit{D}}_A=\mathit{S},$(29) where ${\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}={\mathit{P}}^T{\mathit{X}}_{\mathit{S}}\mathit{P}$ and $\mathit{S}={\mathit{P}}^T({\mathit{A}}^T\mathit{B}+{\mathit{B}}^T\mathit{A})\mathit{P}$. By utilizing the diagonal structure of ${\mathit{D}}_A$, it follows that $x_{i,j}=\frac{s_{ij}}{{\mu }_i+{\mu }_j},$where ${\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}=[x_{i,j}{]}_{i,j=1}^n$. The solution is then given by ${\mathit{X}}_{\mathit{S}}=\mathit{P}{\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}{\mathit{P}}^T.$A necessary and sufficient condition for the uniqueness of ${\mathit{X}}_{\mathit{S}}$, since $\mathit{S}$ is positive definite, is that all the eigenvalues of ${\mathit{A}}^T\mathit{A}$ have a negative real part, that is, ${\mathit{A}}^T\mathit{A}$ is a stable matrix [2,34,35].

4.1. Recursive solution to the symmetric procrustes problem

Consider now the case where the matrices $\mathit{A}$ and $\mathit{B}$ are augmented by the addition of an extra row as in (7(7) $\stackrel{\mathbf{~}}{\mathit{A}}=\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\stackrel{\mathbf{~}}{\mathit{B}}=\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right),$(7) ). The updated symmetric Procrustes problem requires the solution of the optimization problem $\underset{{\stackrel{\mathbf{~}}{\mathit{X}}}_S}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert {\stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}}_S-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}}_S={\stackrel{\mathbf{~}}{\mathit{X}}}_S^T,$where $\stackrel{\mathbf{~}}{\mathit{A}}$ and $\stackrel{\mathbf{~}}{\mathit{B}}$ are defined as in (7(7) $\stackrel{\mathbf{~}}{\mathit{A}}=\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\stackrel{\mathbf{~}}{\mathit{B}}=\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right),$(7) ). Using the Lagrange multipliers and the analysis for the symmetric Procrustes problem as in the previous section, the solution is obtained from the Lyapunov equation (30) ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{X}}}_S+{\stackrel{\mathbf{~}}{\mathit{X}}}_S{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}={\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}+{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$(30) by computing the EVD of ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$. To recursively solve (30(30) ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{X}}}_S+{\stackrel{\mathbf{~}}{\mathit{X}}}_S{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}={\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}+{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$(30) ) observe that on using (28(28) ${\mathit{A}}^T\mathit{A}=\mathit{P}{\mathit{D}}_A{\mathit{P}}^T,$(28) ) yields $\begin{array}{rl}{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}& ={\mathit{A}}^T\mathit{A}+{\mathit{a}}^T\mathit{a}\\ & =\mathit{P}{\mathit{D}}_A{\mathit{P}}^T+{\mathit{a}}^T\mathit{a}\\ & =\mathit{P}({\mathit{D}}_A+{\stackrel{\mathbf{~}}{\mathit{a}}}^T\stackrel{\mathbf{~}}{\mathit{a}}){\mathit{P}}^T.\end{array}$Therefore, the sequential updating of ${\mathit{D}}_A$ requires one rank-1 update, namely, ${\mathit{D}}_A+{\stackrel{\mathbf{~}}{\mathit{a}}}^T\stackrel{\mathbf{~}}{\mathit{a}}={\mathit{P}}_1{\stackrel{\mathbf{~}}{\mathit{D}}}_A{\mathit{P}}_1^T$, whence ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}=\stackrel{\mathbf{~}}{\mathit{P}}{\stackrel{\mathbf{~}}{\mathit{D}}}_A{\stackrel{\mathbf{~}}{\mathit{P}}}^T$with $\stackrel{\mathbf{~}}{\mathit{P}}=\mathit{P}{\mathit{P}}_1$. The Lyapunov equation (30(30) ${\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}{\stackrel{\mathbf{~}}{\mathit{X}}}_S+{\stackrel{\mathbf{~}}{\mathit{X}}}_S{\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{A}}={\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}+{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}}$(30) ) then becomes (31) ${\stackrel{\mathbf{~}}{\mathit{D}}}_A{\stackrel{\mathbf{~}}{\mathit{X}}}_{\mathit{S}}^{\mathit{P}}+{\stackrel{\mathbf{~}}{\mathit{X}}}_{\mathit{S}}^{\mathit{P}}{\stackrel{\mathbf{~}}{\mathit{D}}}_A=\stackrel{\mathbf{~}}{\mathit{S}},$(31) where ${\stackrel{\mathbf{~}}{\mathit{X}}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}={\stackrel{\mathbf{~}}{\mathit{P}}}^T{\stackrel{\mathbf{ˇ}}{\mathit{X}}}_S\stackrel{\mathbf{~}}{\mathit{P}}$ and $\stackrel{\mathbf{~}}{\mathit{S}}={\stackrel{\mathbf{~}}{\mathit{P}}}^T({\stackrel{\mathbf{~}}{\mathit{A}}}^T\stackrel{\mathbf{~}}{\mathit{B}}+{\stackrel{\mathbf{~}}{\mathit{B}}}^T\stackrel{\mathbf{~}}{\mathit{A}})\stackrel{\mathbf{~}}{\mathit{P}}$. The solution to (31(31) ${\stackrel{\mathbf{~}}{\mathit{D}}}_A{\stackrel{\mathbf{~}}{\mathit{X}}}_{\mathit{S}}^{\mathit{P}}+{\stackrel{\mathbf{~}}{\mathit{X}}}_{\mathit{S}}^{\mathit{P}}{\stackrel{\mathbf{~}}{\mathit{D}}}_A=\stackrel{\mathbf{~}}{\mathit{S}},$(31) ) is obtained in a way similar to that for (29(29) ${\mathit{D}}_A{\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}+{\mathit{X}}_{\mathit{S}}^{\mathbf{\left(}\mathit{P}\mathbf{\right)}}{\mathit{D}}_A=\mathit{S},$(29) ).

5. Computational remarks

To amplify the practical usability of the proposed method, the theoretical complexity analysis of the new algorithms has been studied. Consider the general Procrustes problem when no constraint is imposed to the least squares solution matrix $\mathit{X}$. Let the matrices $\stackrel{\mathbf{~}}{\mathit{A}},\stackrel{\mathbf{~}}{\mathit{B}}\in {\mathbb{R}}^{(m+1)\times n}$ as in (7(7) $\stackrel{\mathbf{~}}{\mathit{A}}=\left(\begin{array}{c}\mathit{A}\\ \mathit{a}\end{array}\right)\stackrel{\mathbf{~}}{\mathit{B}}=\left(\begin{array}{c}\mathit{B}\\ \mathit{b}\end{array}\right),$(7) ) and assume that the QRD of $\mathit{A}$ in (6(6) ${\mathit{Q}}_A^T\left(\begin{array}{cc}\mathit{A}& \mathit{B}\end{array}\right)=\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathbf{0}& {\mathit{R}}_{B2}\end{array}\right),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Q}}_A=\left(\begin{array}{cc}{\mathit{Q}}_{A1}& {\mathit{Q}}_{A2}\end{array}\right).$(6) ) is available. The complexity of computing the QRD of $\stackrel{\mathbf{~}}{\mathit{A}}$ afresh is $2n^2(m+1-n/3)$ floating point operations (flops) and that of applying these orthogonal transformations to matrix $\stackrel{\mathbf{~}}{\mathit{B}}$ is $2n^2(2m-n+3)$ flops. The complexities of computing the updating QRD in (9(9) ${\mathit{Q}}_{A_u}^T\left(\begin{array}{cc}{\mathit{R}}_A& {\mathit{R}}_{B1}\\ \mathit{a}& \mathit{b}\end{array}\right)=\left(\begin{array}{cc}{\tilde{\mathit{R}}}_A& {\tilde{\mathit{R}}}_{B_1}\\ \mathbf{0}& {\tilde{\mathit{R}}}_{B_2}\end{array}\right).$(9) ) in order to update ${\mathit{R}}_A$ and ${\mathit{R}}_B$ require $4n^2$ and $8n^2$ flops, respectively. As a result, to compute the QRD of $\stackrel{\mathbf{~}}{\mathit{A}}$ and to apply the orthogonal transformations to $\stackrel{\mathbf{~}}{\mathit{B}}$ will always be computationally more demanding than computing it recursively utilizing previous calculations. The computational efficiency of the recursive method compared to computing the QRD from scratch is, approximately, given by $(9m+2n)/18$.

The solution of the orthogonal Procrustes problem is derived from the EVD of the square matrices $\mathit{F}={\mathit{A}}^T\mathit{B}({\mathit{A}}^T\mathit{B}{)}^T$ and $\mathit{G}={\mathit{B}}^T\mathit{A}{\mathit{A}}^T\mathit{B}$. However, in practice the SVD of ${\mathit{A}}^T\mathit{B}=\mathit{U}\mathbf{\Sigma }{\mathit{V}}^T$ is computed since $\mathit{F}=\mathit{U}{\mathbf{\Sigma }}^2{\mathit{U}}^T=\mathit{U}\mathit{D}{\mathit{U}}^T$ and $\mathit{G}=\mathit{V}{\mathbf{\Sigma }}^2{\mathit{V}}^T=\mathit{V}\mathit{D}{\mathit{V}}^T$. When $\mathit{A}$ and $\mathit{B}$ are modified, the computationally efficient solution to the orthogonal Procrustes problem requires that the SVD of ${\mathit{A}}^T\mathit{B}$ is computed recursively. Therefore, in the recursive solution of the orthogonal and symmetric Procrustes problems, updating SVD decompositions are computed which provide equivalent results as if the EVD were computed. The SVD of an $m\times n$ matrix requires $n^3+m{n}^2+O\left(mn\right)$ flops in the first stage and $km{n}^2+k{n}^2+O\left(mn\right)$ flops in the second stage where k is the number of iterations required to compute each singular value. Herein, it is assumed that $m\ll n$ and it holds that $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{A}}^T\mathit{B}\right)=r<\mathrm{m}\mathrm{i}\mathrm{n}(m,n)$.

The algorithms employ a low-rank modification strategy for the recursive updating of the SVD as in [3]. Algorithms 5.1–5.3 summarize the main computational steps for the solution of the orthogonal Procrustes problem and the recursive proposed solution when new rows or columns of data are added, respectively.

To evaluate further the new algorithms, experiments based on synthetic and real data have been conducted. The computational times of the new algorithms have been compared with the algorithm that solves the same problem afresh in order to obtain the efficiency ratio.

Table 1 presents the execution times (in CPU seconds) of Algorithm 5.1 compared with Algorithm 5.3 when m = 50, 100, 250, 500 rows and n = 1000, 5000 and a single column of data is added to $\mathit{A}$ and $\mathit{B}$, that is k = 1. The results show that keeping n fixed and increasing m reduces the computational efficiency. However, comparing Panel A and B shows that the efficiency is more significant when n increases. Table 2 presents the execution times (in CPU seconds) of the Algorithms 5.1 and 5.2 when m = 100, n = 1000 and when m = 500, n = 10, 000 with a variable number k of new columns are added to $\mathit{A}$ and $\mathit{B}$. All the times presented are the average times after solving the same OPP (afresh or recursively) 100 times. The computational results in Table 2 show that as the number of dimensions increase the computational efficiency also increases. Furthermore, the results suggest that the proposed algorithm for the addition on new column is more efficient when a small number of extra columns is included. Comparing the results in Panel A and B of Table 2 we see that the efficiency of the proposed recursive method increases when the dimensions of the matrices in the original OPP also increase.

Table 1. Execution times in seconds of the recursive solution of the OPP when one single column is added.

	Alg. 5.1	Alg. 5.2	${\displaystyle \frac{\mathrm{A}\mathrm{l}\mathrm{g}.5.1}{\mathrm{A}\mathrm{l}\mathrm{g}.5.2}}$
m	afresh	recursive
(a) Panel A: Orthogonal Procrustes updating with n = 1000 and k = 1.
50	239	6	37
100	273	11	26
250	285	33	9
500	387	178	2
(b) Panel B: Orthogonal Procrustes updating with n = 5000, k = 1.
50	25723	7	3701
100	25510	14	1800
250	25985	71	366
500	26973	307	88

Notes: Results are presented where the matrices in the original Procrustes have m = 50, 100, 250, 500 rows and n = 1000 or 5000 columns and k = 1 column is added.

Table 2. Execution times in seconds of the recursive solution of the OPP when new columns are added.

	Alg. 5.1	Alg. 5.2	${\displaystyle \frac{\mathrm{A}\mathrm{l}\mathrm{g}.5.1}{\mathrm{A}\mathrm{l}\mathrm{g}.5.2}}$
k	afresh	recursive
(a) Panel A: Orthogonal Procrustes updating with m = 100 and n = 1000.
5	2.0	0.1	24
25	2.1	0.1	22
100	2.7	0.2	15
250	3.7	0.5	8
500	6.5	1.6	4
(b) Panel B: Orthogonal Procrustes updating with m = 500 and n = 10, 000.
5	2216.9	11.3	196
25	2117.3	11.3	187
100	2069.8	14.1	147
250	3.7	0.5	8
500	6.5	1.6	4

Notes: Results are presented where the matrices in the original OPP have m = 100 or 500 and n = 1000 and k = 5, 25, 100, 250, 500 new columns are added.

Many a time, the solution of an OPP is applied in the development of algorithms for face recognition, see for example [6,30,37,43]. The algorithm proposed in [43] for feature extraction requires to solve -until convergence- a series of updating OPPs as in (12(12) $\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \mathit{A}\mathit{X}-\mathit{B}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\mathit{X}\mathit{X}}^T={\mathit{X}}^T\mathit{X}=\mathit{I},$(12) ) where $\mathit{A}$ is the data matrix and $\mathit{B}$ is the class indicator matrix after they have both been centred. The processing of an image starts by considering their pixels as the elements of a matrix and then converting each matrix (or each image) to a column vector. Each row in $\mathit{A}$ has length equal to the feature number of the images. When a new image becomes available and is added to the database, a row will be added to $\mathit{A}$. The face recognition algorithm will then have to process the extra image, that is the additional row of the data matrix $\mathit{A}$. This is equivalent to row updating as in (18(18) $\underset{\stackrel{\mathbf{~}}{\mathit{X}}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\Vert \stackrel{\mathbf{~}}{\mathit{A}}\stackrel{\mathbf{~}}{\mathit{X}}-\stackrel{\mathbf{~}}{\mathit{B}}\Vert }_F^2\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{\stackrel{\mathbf{~}}{\mathit{X}}\stackrel{\mathbf{~}}{\mathit{X}}}^T={\stackrel{\mathbf{~}}{\mathit{X}}}^T\stackrel{\mathbf{~}}{\mathit{X}}=\mathit{I},$(18) ). Herein, Algorithm 5.3 is employed to show the computational efficiency of utilizing previous computations when a series of updating OPPs is solved after augmenting the data matrix.

In Table 3 the algorithm that estimates the problems afresh and the recursive algorithm are compared. It is assumed that there are m = 2000, 5000, 10, 000, 15, 000, 20, 000 images available which have been cropped and re-sized to $16\times 16$ and $32\times 32$ pixels, that is n = 256, 1024. The run times (in CPU seconds) when an extra image or alternatively an extra row of data is included 100 times are presented. The results show that when new rows of data are added to the problem the efficiency increases as the number of rows increase. The reason for this is the fact that the recursive algorithm does not depend on m, the number of rows of matrices $\mathit{A}$ and $\mathit{B}$. Instead, it depends on the number of columns n of $\mathit{A}$ and $\mathit{B}$ and on the number of rows l added to the model. This is also why the timings of Algorithm 5.3 are the same when n and l are fixed.

Table 3. Total execution times in seconds of the recursive solution of the OPP when a new single row of data is added 100 times.

	Alg. 5.1	Alg. 5.3	${\displaystyle \frac{\mathrm{A}\mathrm{l}\mathrm{g}.5.1}{\mathrm{A}\mathrm{l}\mathrm{g}.5.3}}$
m	afresh	recursive
(a) Panel A: Orthogonal Procrustes updating with n = 256 and l = 1.
2000	10	6.1	2
5000	20.9	6.1	3
1000	39.3	6.1	6
15,000	57.1	6.1	9
20,000	74.9	6.1	12
(b) Panel B: Orthogonal Procrustes updating with n = 1024 and l = 1.
2000	369.8	371.9	1
5000	546.3	371.9	1
10,000	823.4	371.9	2
15,000	1090.7	371.9	3
20,000	1369.8	371.9	4

Notes: Results are presented when the matrices in the original OPP have dimensions m = 2000, 5000, 10, 000, 15, 000, 20, 000 and n = 256, 1024, and during the updating one extra row is added (l = 1).

Overall, the results show that the computational efficiency of the proposed recursive algorithm increases when the dimensions of the matrices (data) increase. The computational efficiency is significantly more important when a small number of rows or columns is amended -compared to the original matrices- that is, the modification is low-rank. This demonstrates the practicability of the proposed method when solving sequentially OPPs with small modifications in the underlying dataset. The proposed methods are particularly usable when the matrices involved are large-scale and the data are high-dimensional. All the reported computational results were performed on a 64-bit 1.80 GHz Core(TM) i7-8550U processor and 16.0 GB RAM using R (version 3.6.1).

6. Conclusions and future work

The recursive least squares solutions to the matrix equation AX = B when new rows or columns of data are added to data matrices A and B is thoroughly investigated. A computationally efficient algorithm which is based on the singular value decomposition is proposed. Computational results are presented for synthetic data and also for a machine learning application based on feature extraction for face recognition. The recursive solution to the symmetric Procrustes problem when including extra data is also investigated. The experimental results suggest that the proposed algorithms are computationally more efficient when the matrices are high-dimensional and when they are augmented with a small number of rows or columns.

The extension of the proposed method to include other special classes of matrices like reflexive and anti-reflexive, Stiefel matrices or Toeplitz matrices merits further investigation. Future work will also consider the solution to the Procrustes problem with regularization constraints after modifying the matrices with the addition or deletion of data.

Disclosure statement

No potential conflict of interest was reported by the author(s).