# A theory of meta-factorization

Michał P. Karpowicz
Warsaw University of Technology
and
NASK National Research Institute
<EMAIL_ADDRESS>
###### Abstract

We introduce meta-factorization, a theory that describes matrix decompositions
as solutions of linear matrix equations: the projector and the reconstruction
equation. Meta-factorization reconstructs known factorizations, reveals their
internal structures, and allows for introducing modifications, as illustrated
with SVD, QR, and UTV factorizations. The prospect of meta-factorization also
provides insights into computational aspects of generalized matrix inverses
and randomized linear algebra algorithms. The relations between the Moore-
Penrose pseudoinverse, generalized Nyström method, the CUR decomposition, and
outer-product decomposition are revealed as an illustration. Finally, meta-
factorization offers hints on the structure of new factorizations and provides
the potential of creating them.

_K_ eywords Matrix factorization; Generalized inverses; Low-rank matrix
approximation; Randomized algorithms

##  1 Introduction

Matrix factorization decomposes a matrix into a product of other matrices to
reveal essential properties of the original matrix. We introduce meta-
factorization, a procedure that further factorizes the factors to reveal
details of the original decomposition and control its internal structure. The
main asset of meta-factorization is the projector equation, which reveals
hidden projectors of matrix decomposition and describes the mutual
interactions of the factors needed to reconstruct the original matrix. The
reconstruction equation then governs the reconstruction.

Identifying the hidden projectors is one of the key benefits of meta-
factorization. Once their internal structure is revealed, the projector
equation keeps the reconstruction abilities of the projectors intact while
allowing for modifications to be introduced when necessary. The original
matrix is safely mapped onto its column and row space. As a result, meta-
factorization develops a nontrivial perspective in which matrix factorizations
emerge as solutions to linear matrix equations.

This prospect becomes useful when studying matrix decompositions and designing
matrix algorithms, including randomized linear algebra and low-rank
approximation algorithms. By following the rules of meta-factorization, we can
replace selected factors of a given decomposition with their numerically
attractive equivalents, which control the performance and stability of the
design. Finally, meta-factorization provides hints on the structure of new
factorizations and offers the potential to create them.

###  1.1 Related work

Matrix factorizations play a critical role in mathematics, science, and
engineering. They reveal essential matrix properties, such as spectrum, rank,
or fundamental subspaces, which in turn characterize solutions to many
fundamental problems in physics, signal processing, control engineering,
artificial intelligence, and data science, see, e.g., Strang [30, 29]. Matrix
factorizations also correspond to matrix algorithms, which is clearly
illustrated by Golub and Van Loan [8], Stewart [28], Demmel [5], or Trefethen
[32, 33]. Moreover, Edelman and Jeong [6] show that there are still new
families of matrix factorizations to be analyzed together with their unifying
algebraic structure. The search for new matrix factorizations has been fueling
linear algebra for years.

From a mathematical point of view, meta-factorization is related to the theory
of generalized matrix inverses. The equations of meta-factorization lead to
solutions involving pseudoinverses or define them in special cases. For the
related fundamental results, see Penrose [24] and Langenhop [15], Matsaglia
and Styan [21] and Ben-Israel and Greville [1] for the comprehensive
exposition. Meta-factorization refers to computational aspects of generalized
inverses as well. In particular, it describes explicit formulas for the Moore-
Penrose pseudoinverse. Such formulas, critical for the accuracy and stability
of computations, have been summarized by Ben-Israel and Greville in [1]. See
also James [13] for an introduction, and Strang and Drucker [31] for a fresh
perspective. One more relation should be mentioned when a big picture of
matrix factorizations is considered. Namely, there is a connection between the
reconstruction equation of meta-factorization and the quadratic matrix
equation generating the Lie groups. The latter, studied by Edelman and Jeong
in [6], gives rise to families of matrix factorizations. Under the appropriate
assumptions, it can be interpreted as the meta-factorization equation.

The results presented in this paper were initially inspired by the work of
Sorensen and Embree [27], with their technique for deriving the column and row
space projectors for the CUR decomposition. This idea in the form of a
projector equation is included in the process of meta-factorization. Also, as
we will see, the generalized form of the projector equation may pave the way
to possible new factorizations.

Designing projectors has been of great interest in the field of randomized
linear algebra, which is visible in the works of Halko et al. [10], Mahoney
[16], Martinsson and Tropp [18, 20], Shabat et al. [26], and Kannan and
Vempala [14]. Recently, Nakatsukasa [22] showed that many algorithms of
randomized linear algebra could be written in the form of the generalized
Nyström’s method, a two-sided projection defined by sampling (sketch)
matrices. Meta-factorization gives rise to the generalized Nyström’s method.
It also explains the structure of the CUR decomposition studied by Mahoney and
Drineas in [17], Wang and Zhang [34], Sorensen and Embree [27], Hamm and Huang
[11], and predicted by Penrose in [25]. The outer-product decomposition
resulting from the the rank-reduction process introduced by Wedderburn [35],
Egerváry [7], Householder [12], Cline and Funderlic [4], and Chu et al. [3],
is also described by the meta-factorization equations. Finally, Nakatsukasa
and Tropp in [23] show remarkable results of redesigning projectors in GMRES
and Rayleigh-Ritz methods. In particular, they predict the benefits of moving
into non-orthogonal linear algebra, a move meta-factorization fully supports.

###  1.2 Outline

The paper is divided into four sections. Section 2 introduces the concept of
meta-factorization. First, it demonstrates how matrix factors interact to
reconstruct the original matrix. Next, it describes the nature of those
interactions in the form of a projection equation and defines the equations of
meta-factorizations. Finally, it shows that solutions of the equations define
matrix factorizations. Section 3 shows how meta-factorization works. The
theory is applied to reconstruct the singular value decomposition (SVD),
examine the internal structure of column-pivoted QR decomposition (CPQR), and
develop a family of UTV factorizations. Section 4 discusses additional
benefits of meta-factorization. It presents a relationship of meta-
factorization with the theory of generalized matrix inverses, reconstructs the
generalized Nyström’s method for computing low-rank matrix approximations, and
shows its relationship with the CUR decomposition111A Matlab code supplement
is available here: https://github.com/mkarpowi/mft. Finally, further research
directions are proposed in the Appendix A.

#### Notation.

Capital bold letters denote matrices, $\mathbf{W}^{*}$ denotes the conjugate
transpose of $\mathbf{W}$, and $\mathbf{W}^{+}$ denotes the (Moore-Penrose
inverse) pseudoinverse of $\mathbf{W}$. The $k\times k$ identity matrix is
denoted as $\mathbf{I}_{k}$. Given ordered subindex sets $I$ and $J$,
$\mathbf{A}(I,J)$ denotes the submatrix of $\mathbf{A}$ containing rows and
columns of $\mathbf{A}$ indexed by $I$ and $J$, with $\mathbf{A}(:,J)$
extracting the columns of $\mathbf{A}$ indexed by $J$. With $k$ being a
positive integer, $1{:}k$ denotes the ordered set $(1,\dots,k)$. Finally,
$\mathcal{C}(\mathbf{A})$ denotes the column space and
$\mathcal{C}(\mathbf{A}^{*})$ denotes the row space for $\mathbf{A}$. An
orthonormal $m\times n$ matrix $\mathbf{U}$ has orthonormal columns:
$\mathbf{U}^{*}\mathbf{U}=\mathbf{I}_{n}$. A unitary matrix is a square
orthonormal matrix.

##  2 Meta-factorization

This section introduces the theoretical foundations of the proposed concept of
meta-factorization. First, by investigating the internal structure of certain
matrix factorizations, we demonstrate how the factors interact with each other
to reconstruct the original matrix. Next, we derive the conditions describing
the nature of those interactions. Finally, the concept of meta-factorization
is formulated.

###  2.1 Factorizations and hidden projections

The matrix factorizations we study in this paper are of the following form:

$\displaystyle\setstackgap{L}{15pt}$

where an arbitrary $m$ by $n$ matrix $\mathbf{A}$ is decomposed here into the
product of an $m$ by $k$ matrix $\mathbf{F}$, a square $k$ by $k$ matrix
$\mathbf{G}$, and an $n$ by $k$ matrix $\mathbf{H}$. These factors may reveal
essential matrix properties, such as spectrum, rank, or four fundamental
subspaces. The most important factorizations, such as the singular value
decomposition (SVD), the column-pivoted QR (CPQR) decomposition, the
eigenvalue decomposition (ED), the LU decomposition, or the CUR decomposition,
are all described by this pattern.

Matrix $\mathbf{F}$ defines a basis for the column space, whereas $\mathbf{H}$
defines a basis for the row space of $\mathbf{A}$. Those two bases can be
derived directly from the columns and rows of the original matrix, for example
as a result of (random) sampling, elimination, or orthogonalizing
transformations. Matrix $\mathbf{G}$ turns out to be a more challenging one.
It encodes the essential properties of $\mathbf{A}$ isolated by $\mathbf{F}$
and $\mathbf{H}$ in a way that guarantees their reconstruction in the original
matrix.

The key question is how to determine the mixing matrix $\mathbf{G}$ when
$\mathbf{F}$ and $\mathbf{H}$ are given? As mentioned, the best way to learn
more about a given matrix is to factor it into a particular mixture of other
matrices. Let us apply the same trick here: let us factor the mixing matrix
$\mathbf{G}$.

A straightforward way to calculate $\mathbf{G}$ as a product of other matrices
is to invert (2.1):

$\displaystyle\mathbf{G}=\mathbf{F}^{+}\mathbf{A}(\mathbf{H}^{*})^{+}.$ (1)

However, this operation does not reveal everything there is to reveal.
Fortunately, there is a more promising way. Instead of inverting $\mathbf{F}$
and $\mathbf{H}$, we can introduce new input and output transformations,
$\mathbf{Y}$ and $\mathbf{X}$, and assume that

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{G}}{\scriptstyle(k\times
k)}&=\stackunder{\mathbf{Y}^{*}}{\scriptstyle(k\times
m)}&&\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&&\stackunder{\mathbf{X}.}{\scriptstyle(n\times k)}\end{aligned}$ (2)

The key insight follows from substituting (2) into (2.1):

$\displaystyle\mathbf{A}=\mathbf{F}\mathbf{Y}^{*}\mathbf{A}\mathbf{X}\mathbf{H}^{*}.$
(3)

The original decomposition (2.1) is a reconstruction formula in disguise.
Matrix factorization reconstructs the original matrix by projecting it onto
its own column space and row space. For that to be possible it defines two
projectors, the column space projector $\mathbf{F}\mathbf{Y}^{*}$ and the row
space projector $\mathbf{X}\mathbf{H}^{*}$. Those projectors are determined by
$\mathbf{F}$ and $\mathbf{H}$, but also by $\mathbf{Y}$ and $\mathbf{X}$, the
factors encoding the properties of $\mathbf{A}$ in $\mathbf{G}$. Therefore, to
calculate the mixing matrix when $\mathbf{F}$ and $\mathbf{H}$ are known, it
is necessary to find $\mathbf{Y}$ and $\mathbf{X}$ that define projectors for
the column space and the row space of $\mathbf{A}$. That is the idea behind
meta-factorization.

###  2.2 The projector equation

As we have seen, matrix factorization reconstructs the original matrix by
projecting it onto its own column and row space. Matrix decomposition given by
(2.1) is a reconstruction formula (3) in disguise. The idea behind meta-
factorization is to control (or shape) the internal structure of that
reconstruction formula. This, however, can only be accomplished provided that
the properties of the projection matrices remain intact. We formulate that
constraint in the form of a projector equation.

Suppose that $\mathbf{F}$ and $\mathbf{H}$ are arbitrary matrices. The
projector equation defines $\mathbf{Y}$ and $\mathbf{X}$ such that

$\displaystyle\mathbf{P}=\mathbf{F}\mathbf{Y}^{*}\quad\text{and}\quad\mathbf{R}=\mathbf{X}\mathbf{H}^{*}$
(4)

are idempotent matrices, i.e., $\mathbf{P}^{2}=\mathbf{P}$ and
$\mathbf{R}^{2}=\mathbf{R}$.

######  Theorem 1.

If $\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$ are idempotent
rank-$k$ matrices, then they satisfy the projector equation,

$\displaystyle\mathbf{Y}^{*}\mathbf{F}=\mathbf{H}^{*}\mathbf{X}=\mathbf{I}_{k}.$
(5)

Conversely, if the the projector equation holds, then
$\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$ are idempotent
rank-$k$ matrices.

###### Proof.

Consider $\mathbf{P}=\mathbf{F}\mathbf{Y}^{*}$ and
$\mathbf{R}=\mathbf{X}\mathbf{H}^{*}$. If (5) holds, then
$\mathbf{P}^{2}=\mathbf{F}\mathbf{Y}^{*}\mathbf{F}\mathbf{Y}^{*}=\mathbf{F}\mathbf{Y}^{*}=\mathbf{P}$.
By assumption, $\mathrm{rank}(\mathbf{F})\leq k$, so we have

$\displaystyle\mathrm{rank}(\mathbf{F}\mathbf{Y}^{*})\leq\mathrm{rank}(\mathbf{F})\leq
k=\mathrm{rank}(\mathbf{Y}^{*}\mathbf{F})\leq\mathrm{rank}(\mathbf{F})=\mathrm{rank}(\mathbf{F}\mathbf{Y}^{*}\mathbf{F})\leq\mathrm{rank}(\mathbf{F}\mathbf{Y}^{*}),$
(6)

so $\mathrm{rank}(\mathbf{P})=k$. Similarly,
$\mathbf{R}^{2}=\mathbf{X}\mathbf{H}^{*}\mathbf{X}\mathbf{H}^{*}=\mathbf{X}\mathbf{H}^{*}=\mathbf{R}$,
and $\mathrm{rank}(\mathbf{R})=k$, by the same arguments.

Conversely, if $\mathbf{F}\mathbf{Y}^{*}$ is idempotent rank-$k$ matrix, then

$\displaystyle\mathbf{F}(\mathbf{Y}^{*}\mathbf{F}-\mathbf{I}_{k})\mathbf{Y}^{*}=\mathbf{0}$
(7)

and $\mathrm{rank}(\mathbf{F})=\mathrm{rank}(\mathbf{Y}^{*})=k$. Indeed,
suppose $\mathbf{F}$ has linearly dependent columns and the equation holds for
a nonzero matrix $(\mathbf{Y}^{*}\mathbf{F}-\mathbf{I}_{k})\mathbf{Y}^{*}$.
But then $\mathbf{F}\mathbf{Y}^{*}$ is not a rank-$k$ matrix, which
contradicts the assumption. The same argument holds for $\mathbf{Y}^{*}$.
Therefore, we conclude that (5) holds.

Similarly, if $\mathbf{X}\mathbf{H}^{*}$ is idempotent rank-$k$, then

$\displaystyle\mathbf{X}(\mathbf{H}^{*}\mathbf{X}-\mathbf{I}_{k})\mathbf{H}^{*}=\mathbf{0}.$
(8)

Since $\mathrm{rank}(\mathbf{X})=\mathrm{rank}(\mathbf{H}^{*})=k$, we conclude
that (5) holds. ∎

We shall call equation (5) the projector equation. The equation is fundamental
and paves the way to our meta-factorization results. It connects
$\mathbf{Y}^{*}$ with $\mathbf{F}$ and $\mathbf{X}$ with $\mathbf{H}^{*}$ to
make them projectors, which is a key step of the meta-factorization procedure.

The equation has been often used in linear algebra. For example, it appears in
the already mentioned works of Langenhop [15], Ben-Israel and Greville [1] on
the generalized inverses of matrices, and in the recent work of Sorensen and
Embree [27] on the CUR decomposition. However, it seems that its role has not
been extensively explored in the context being considered. The projector
equation can also be generalized, in which case it may give rise to new
factorizations, as predicted in Appendix A.

To appreciate the role of the projector equation, we must characterize its
solutions.

######  Theorem 2.

The projector equation (5) admits solutions:

$\displaystyle\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}$
(9)

for any $\mathbf{B}$ and $\mathbf{D}$ such that:

$\displaystyle\mathrm{rank}(\mathbf{B}^{*}\mathbf{F})=\mathrm{rank}(\mathbf{H}^{*}\mathbf{D})=k=\min\\{m,n\\}.$
(10)

Furthermore, $\mathbf{Y}^{*}$ and $\mathbf{X}$ satisfy the generalized matrix
inverse equations:

$\displaystyle\mathbf{F}\mathbf{Y}^{*}\mathbf{F}=\mathbf{F}\quad\text{and}\quad\mathbf{H}^{*}\mathbf{X}\mathbf{H}^{*}=\mathbf{H}^{*}.$
(11)

Finally, when $\mathbf{B}=\mathbf{F}$ and $\mathbf{D}=\mathbf{H}$, then
$\mathbf{Y}^{*}=\mathbf{F}^{+}$ and $\mathbf{X}=(\mathbf{H}^{*})^{+}$.

###### Proof.

Substituting (9) into (5), we obtain

$\displaystyle\mathbf{Y}^{*}\mathbf{F}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}=\mathbf{I}_{k}=\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}=\mathbf{H}^{*}\mathbf{X}.$
(12)

It follows that $\mathbf{Y}^{*}$ and $\mathbf{X}$ satisfy the generalized
inverse equations (11).

To prove that $\mathbf{Y}^{*}=\mathbf{F}^{+}$ when $\mathbf{B}=\mathbf{F}$, we
show that in the general case,

$\displaystyle(\mathbf{F}^{*}\mathbf{F})^{+}\mathbf{F}^{*}=\mathbf{F}^{+}.$
(13)

Let $\mathbf{W}=(\mathbf{F}^{*}\mathbf{F})^{+}\mathbf{F}^{*}$. We must verify
that $\mathbf{W}$ satisfies the Penrose equations [24] for a pseudoinverse:

$\displaystyle\begin{aligned} \mathbf{F}\mathbf{W}\mathbf{F}&=\mathbf{F},\\\
\mathbf{W}\mathbf{F}\mathbf{W}&=\mathbf{W},\\\
(\mathbf{F}\mathbf{W})^{*}&=\mathbf{F}\mathbf{W},\\\
(\mathbf{W}\mathbf{F})^{*}&=\mathbf{W}\mathbf{F}.\end{aligned}$ (14)

First, consider $\mathbf{F}^{*}\mathbf{F}$ and observe that by the properties
of a pseudoinverse,

$\displaystyle\mathbf{F}^{*}\mathbf{F}=\mathbf{F}^{*}\mathbf{F}(\mathbf{F}^{*}\mathbf{F})^{+}\mathbf{F}^{*}\mathbf{F}=\mathbf{F}^{*}\mathbf{F}\mathbf{W}\mathbf{F},$
(15)

which implies that

$\displaystyle\mathbf{F}^{*}\mathbf{F}(\mathbf{W}\mathbf{F}-\mathbf{I}_{k})=\mathbf{0}.$
(16)

Multiplying on the left by $(\mathbf{W}\mathbf{F}-\mathbf{I}_{k})^{*}$ and
using the properties of transposition gives
$(\mathbf{F}\mathbf{W}\mathbf{F}-\mathbf{F})^{*}(\mathbf{F}\mathbf{W}\mathbf{F}-\mathbf{F})=\mathbf{0}$.
Therefore,

$\displaystyle\mathbf{F}\mathbf{W}\mathbf{F}=\mathbf{F}$ (17)

and the first Penrose equation is satisfied.

Next, consider the equation,

$\displaystyle(\mathbf{F}^{*}\mathbf{F})^{+}\mathbf{F}^{*}\mathbf{F}(\mathbf{F}^{*}\mathbf{F})^{+}=(\mathbf{F}^{*}\mathbf{F})^{+}.$
(18)

Multiplying on the right by $\mathbf{F}^{*}$, we conclude that

$\displaystyle\mathbf{W}\mathbf{F}\mathbf{W}=\mathbf{W}.$ (19)

Therefore, the second Penrose equation is satisfied, as well.

The last two Penrose equations follow from the properties of transpositions.
Namely, we have

$\displaystyle\mathbf{F}\mathbf{W}=(\mathbf{F}\mathbf{W})^{*}\quad\text{and}\quad\mathbf{W}\mathbf{F}=(\mathbf{W}\mathbf{F})^{*}$
(20)

and we conclude that

$\displaystyle\mathbf{W}=\mathbf{F}^{+}.$ (21)

By similar arguments,

$\displaystyle(\mathbf{H}^{*})^{+}=\mathbf{H}(\mathbf{H}^{*}\mathbf{H})^{+}$
(22)

and the result holds with $\mathbf{D}=\mathbf{H}$. ∎

We have demonstrated that matrices of the form given by (9) solve the
projector equation and give rise to orthogonal projectors when
$\mathbf{B}=\mathbf{F}$ and $\mathbf{D}=\mathbf{H}$. It remains to show that
they also define oblique projectors. In other words, $\mathbf{Y}^{*}$ and
$\mathbf{X}$ must satisfy the generalized inverse equations even when
$\mathbf{B}^{*}\mathbf{F}$ and $\mathbf{H}^{*}\mathbf{D}$ are not full rank
matrices.

######  Theorem 3.

Suppose that

$\displaystyle\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+},$
(23)

where $\mathbf{B}$ and $\mathbf{D}$ satisfy condition:

$\displaystyle\mathrm{rank}(\mathbf{B}^{*}\mathbf{F})=\mathrm{rank}(\mathbf{H}^{*}\mathbf{D})=k.$
(24)

Then

$\displaystyle\mathbf{F}\mathbf{Y}^{*}\mathbf{F}=\mathbf{F}\quad\text{and}\quad\mathbf{H}^{*}\mathbf{X}\mathbf{H}^{*}=\mathbf{H}^{*}.$
(25)

###### Proof.

It follows from the properties of pseudoinverse that:

$\displaystyle\mathbf{B}^{*}\mathbf{F}(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}-\mathbf{B}^{*}\mathbf{F}=\mathbf{0}.$
(26)

Multiplying on the left by
$((\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}-\mathbf{I}_{k})^{*}$,
yields

$\displaystyle(\mathbf{B}(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}-\mathbf{B})^{*}(\mathbf{F}(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}-\mathbf{F})=\mathbf{0}.$
(27)

Notice that $\mathbf{B}^{*}\mathbf{F}$ is not a full rank matrix, so we cannot
expect $(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}$ to be equal to
$\mathbf{I}_{k}$. Therefore, for the equation to hold for every $\mathbf{B}$,
we must have

$\displaystyle\mathbf{F}(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{F}=\mathbf{F}\mathbf{Y}^{*}\mathbf{F}=\mathbf{F}.$
(28)

Similarly, consider the equation

$\displaystyle\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{H}^{*}\mathbf{D}-\mathbf{H}^{*}\mathbf{D}=(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}-\mathbf{I}_{k})\mathbf{H}^{*}\mathbf{D}=\mathbf{0}.$
(29)

Multiplying on the right by
$(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}-\mathbf{I}_{k})^{*}$,
yields

$\displaystyle\begin{aligned}
&(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}-\mathbf{I}_{k})\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}-\mathbf{I}_{k})^{*}\\\
&=(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{H}^{*}-\mathbf{H}^{*})\mathbf{D}^{**}(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}-\mathbf{I}_{k})^{*}\\\
&=(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{H}^{*}-\mathbf{H}^{*})(\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{D}^{*}-\mathbf{D}^{*})^{*}=\mathbf{0}.\end{aligned}$
(30)

Again, since $\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}$ is not
equal to $\mathbf{I}_{k}$ in the general case considered and the equation must
hold for every $\mathbf{D}$, we conclude that

$\displaystyle\mathbf{H}^{*}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{H}^{*}=\mathbf{H}^{*}\mathbf{X}\mathbf{H}^{*}=\mathbf{H}^{*}.$
(31)

∎

In short, the projector equation defines necessary and sufficient conditions
for $\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$ to be idempotent
rank-$k$ matrices. The form of the solutions $\mathbf{Y}^{*}$ and
$\mathbf{X}$, given by (9), can be generalized to make
$\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$ oblique projectors
(idempotent matrices).

###  2.3 The reconstruction equation and meta-factorization procedure

We are now ready to introduce the concept of meta-factorization, a procedure
for decomposing the factorization’s factors into the product of matrices of
the desired structure.

When $\mathbf{F}$ and $\mathbf{H}$ are given, the projector equation defines
matrices $\mathbf{Y}$ and $\mathbf{X}$ that satisfy the constraints of the
matrix reconstruction problem. It keeps the properties of the projection
matrices, $\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$, intact.
Therefore, we must also enforce the projections onto the required subspaces.
Indeed, notice that if the projector equation holds, $\mathbf{A}$ is
transformed into some

$\displaystyle\mathbf{M}=(\mathbf{F}\mathbf{Y}^{*})\mathbf{A}(\mathbf{X}\mathbf{H}^{*}).$
(32)

But when $\mathbf{F}$ and $\mathbf{H}$ are arbitrary, we cannot expect
$\mathbf{M}$ to approximate $\mathbf{A}$. For the reconstruction to be
successful it is also necessary to demand that $\mathbf{F}$ and $\mathbf{H}$
define the column space and row space of $\mathbf{A}$. That brings us to our
main result. Theorem 4 introduces the idea of meta-factorization exploiting
the projector equation to find and tune the desired projector matrices.

######  Theorem 4.

The reconstruction equation

$\displaystyle\mathbf{A}=\mathbf{F}\mathbf{Y}^{*}\mathbf{A}\mathbf{X}\mathbf{H}^{*}$
(33)

holds for a rank-$k$ matrix $\mathbf{A}$, when $\mathbf{Y}^{*}$ and
$\mathbf{X}$ solve the projector equation (5) with $\mathbf{F}$ and
$\mathbf{H}$ representing bases for the column and row space of $\mathbf{A}$.

###### Proof.

The result follows from the elementary properties of projections. By Theorem
1, $\mathbf{P}=\mathbf{F}\mathbf{Y}^{*}$ is a projection. By the assumption
made, we also have $\mathcal{C}(\mathbf{A})=\mathcal{C}(\mathbf{P})$. Take any
$\mathbf{a}\in\mathcal{C}(\mathbf{A})=\\{\mathbf{y}\in\mathbb{R}^{m}\mid\mathbf{y}=\mathbf{P}\mathbf{x}\\}$.
It follows that

$\displaystyle\mathbf{P}\mathbf{a}=\mathbf{P}(\mathbf{P}\mathbf{x})=\mathbf{P}\mathbf{x}=\mathbf{a}.$
(34)

The same argument holds for $\mathbf{R}=\mathbf{X}\mathbf{H}^{*}$. For any
$\mathbf{b}\in\mathcal{C}(\mathbf{A}^{*})=\\{\mathbf{w}\in\mathbb{R}^{n}\mid\mathbf{w}=\mathbf{R}^{*}\mathbf{v}\\}$,
we have

$\displaystyle\mathbf{R}^{*}\mathbf{b}=\mathbf{b}.$ (35)

As a result, $\mathbf{P}\mathbf{A}=\mathbf{A}$ and
$\mathbf{R}^{*}\mathbf{A}^{*}=\mathbf{A}^{*}$, which gives
$\mathbf{P}\mathbf{A}\mathbf{R}=\mathbf{A}$. ∎

Basic meta-factorization proceeds according to the following steps
(illustrated with Matlab code):

  1. Step 1:

Given a matrix $\mathbf{A}$, select basis $\mathbf{F}$ for the column space
$\mathcal{C}(\mathbf{A})$ and basis $\mathbf{H}$ for the row space
$\mathcal{C}(\mathbf{A}^{*})$.

  2. Step 2:

Solve the projector equation
$\mathbf{Y}^{*}\mathbf{F}=\mathbf{H}^{*}\mathbf{X}=\mathbf{I}_{k}$ to find the
input basis $\mathbf{X}$ and the output basis $\mathbf{Y}$:

        [QY,RY] = qr(B’*F,0);       % If needed, select B and D,
    [QX,RX] = qr(H’*D,0);       % otherwise, set B = F and D = H.
    Y       = (RY\(QY’*B’))’;   % If rank(B’*F) = rank(H’*D) = k,
    X       = (D/RX)*QX’;       % calculate Y and X.

  3. Step 3:

Calculate the mixing matrix $\mathbf{G}=\mathbf{Y}^{*}\mathbf{A}\mathbf{X}$
and (if needed) perform its factorization to obtain the desired decomposition
$\mathbf{A}=\mathbf{F}\mathbf{G}\mathbf{H}^{*}$:

        G       = Y’*A*X;
    % If needed, perform additional factorization of G, e.g.,
    % [U,S,V] = svd(G), [Q,R] = qr(G) or [L,U] = lu(G), etc.
    Ar      = F*G*H’;

The procedure finds solutions of the meta-factorization equations,

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{G}}{\scriptstyle(k\times
k)}&=\stackunder{\mathbf{Y}^{*}}{\scriptstyle(k\times
m)}&&\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&&\stackunder{\mathbf{X},}{\scriptstyle(n\times k)}\\\
\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&=\stackunder{\mathbf{F}}{\scriptstyle(m\times
k)}&&\stackunder{\mathbf{G}}{\scriptstyle(k\times
k)}&&\stackunder{\mathbf{H}^{*}.}{\scriptstyle(k\times n)}\end{aligned}$ (36)

Matrices $\mathbf{F}$ and $\mathbf{H}$ are the parameters of the procedure,
and their role is to define bases for the column space and the row space of
$\mathbf{A}$. Given $\mathbf{F}$ and $\mathbf{H}$, matrices $\mathbf{X}$,
$\mathbf{Y}$, and $\mathbf{G}$ are determined in a way that allows
$\mathbf{A}$ to be reconstructed from $\mathbf{G}$. As a result, $\mathbf{A}$
is decomposed into the product of $\mathbf{F}$, $\mathbf{G}$, and
$\mathbf{H}$, and $\mathbf{G}$ into the product of $\mathbf{Y}$, $\mathbf{X}$
and $\mathbf{A}$. Finally, factorization of the mixing matrix $\mathbf{G}$,
i.e., meta-factorization, shapes the structure of projectors and their
numerical properties.

Notice that the proof of Theorem 4 only requires idempotent
$\mathbf{F}\mathbf{Y}^{*}$ and $\mathbf{X}\mathbf{H}^{*}$. We can take
advantage of that fact and use Theorem 3 to formulate a more general result.

######  Corollary 1.

The reconstruction equation holds, if

$\displaystyle\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+},$
(37)

where $\mathbf{B}$ and $\mathbf{D}$ satisfy condition:

$\displaystyle\mathrm{rank}(\mathbf{B}^{*}\mathbf{F})=\mathrm{rank}(\mathbf{H}^{*}\mathbf{D})=k,$
(38)

and $\mathbf{F}$ and $\mathbf{H}$ represent bases for the column and row space
of $\mathbf{A}$.

We can further adjust the meta-factorization procedure to admit low-rank bases
as well. Given a rank-$k$ matrix $\mathbf{A}$ and some positive $r\leq k$, we
can find the required rank-$r$ projectors by solving the projector equation
(5) with $\mathbf{I}_{r}$ for the low-rank bases $\mathbf{F}_{r}$ and
$\mathbf{H}_{r}$. Setting $\mathbf{B}=\mathbf{F}_{r}$ and
$\mathbf{D}=\mathbf{H}_{r}$ is preferred in this case. The procedure also
admits designing one projector only, in which case it is enough to work with
either $\mathbf{F}$ or $\mathbf{H}$.

###  2.4 Factorizations as solutions of linear matrix equations

Finally, this section develops the prospect of meta-factorization that
describes matrix factorizations as solutions of linear matrix equations.

In his seminal work on the general inverse for matrices, Penrose characterized
solutions to linear matrix equations. The statement of the theorem is
presented below in the notation of the present paper.

######  Theorem 5 (Penrose [24]).

A necessary and sufficient condition for the equation

$\displaystyle\mathbf{F}\mathbf{G}\mathbf{H}^{*}=\mathbf{A}$ (39)

to have a solution $\mathbf{G}$ is

$\displaystyle\mathbf{F}\mathbf{F}^{+}\mathbf{A}(\mathbf{H}^{*})^{+}\mathbf{H}^{*}=\mathbf{A},$
(40)

in which case the general solution is

$\displaystyle\mathbf{G}=\mathbf{F}^{+}\mathbf{A}(\mathbf{H}^{*})^{+}+\mathbf{W}-\mathbf{F}^{+}\mathbf{F}\mathbf{W}\mathbf{H}^{*}(\mathbf{H}^{*})^{+},$
(41)

where $\mathbf{W}$ is arbitrary.

The proof presented by Penrose exploits the properties of pseudoinverses,
introduced in the very same paper. Meta-factorization, exploiting the
properties of projections, shows that another perspective can be taken in
which solutions to the addressed matrix equation become matrix factorizations.
Indeed, when we interpret (39) as a description of a matrix factorization,
then equation (40) defines that matrix factorization by the properly
constructed projection matrices. The latter equation is clearly a special case
of the reconstruction equation (33) of Theorem 4. In other words, matrix
factorization given by (39) is a solution of the meta-factorization equations
(33).

######  Theorem 6.

Given $\mathbf{A}$, $\mathbf{F}$ and $\mathbf{H}$, a necessary and sufficient
condition for the equation

$\displaystyle\mathbf{F}\mathbf{G}\mathbf{H}^{*}=\mathbf{A}$ (42)

to have a solution $\mathbf{G}$ is

$\displaystyle\mathbf{F}\mathbf{Y}^{*}\mathbf{A}\mathbf{X}\mathbf{H}^{*}=\mathbf{A},$
(43)

where

$\displaystyle\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+},$
(44)

and
$\mathrm{rank}(\mathbf{B}^{*}\mathbf{F})=k=\mathrm{rank}(\mathbf{H}^{*}\mathbf{D})$.
The general solution is

$\displaystyle\mathbf{G}=\mathbf{Y}^{*}\mathbf{A}\mathbf{X}+\mathbf{W}-\mathbf{Y}^{*}\mathbf{F}\mathbf{W}\mathbf{H}^{*}\mathbf{X},$
(45)

where $\mathbf{W}$ is arbitrary.

###### Proof.

If $\mathbf{G}$ satisfies (42), then
$\mathcal{C}(\mathbf{F})=\mathcal{C}(\mathbf{A})$ and
$\mathcal{C}(\mathbf{H}^{*})=\mathcal{C}(\mathbf{A}^{*})$. By Theorem 3 and
Corollary 1 we have

$\displaystyle\mathbf{A}$
$\displaystyle=\mathbf{F}(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{A}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}\mathbf{H}^{*},\quad\text{where}\quad\mathbf{G}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{A}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}=\mathbf{Y}^{*}\mathbf{A}\mathbf{X}.$
(46)

Conversely, if (43) holds, then
$\mathbf{G}=(\mathbf{B}^{*}\mathbf{F})^{+}\mathbf{B}^{*}\mathbf{A}\mathbf{D}(\mathbf{H}^{*}\mathbf{D})^{+}$
is a particular solution of (42). The general solution must also satisfy the
homogeneous equation, $\mathbf{F}\mathbf{G}\mathbf{H}^{*}=\mathbf{0}$.
Consider

$\displaystyle\mathbf{G}=\mathbf{W}-\mathbf{Y}^{*}\mathbf{F}\mathbf{W}\mathbf{H}^{*}\mathbf{X}.$
(47)

By Theorem 2 and 3, we have

$\displaystyle\begin{aligned}
\mathbf{F}\mathbf{G}\mathbf{H}^{*}&=\mathbf{F}\mathbf{W}\mathbf{H}^{*}-\mathbf{F}\mathbf{Y}^{*}\mathbf{F}\mathbf{W}\mathbf{H}^{*}\mathbf{X}\mathbf{H}^{*}\\\
&=\mathbf{F}\mathbf{W}\mathbf{H}^{*}-\mathbf{F}\mathbf{W}\mathbf{H}^{*}=\mathbf{0}.\end{aligned}$
(48)

Therefore, the general solution to (42) is given by (45). ∎

Notice that when the projector equation holds, the solution of the homogeneous
equation is $\mathbf{G}=\mathbf{0}$. Otherwise, when the projector equation is
not satisfied, $\mathbf{G}$ is rectangular and depends on $\mathbf{W}$.

As a special case of the result above we get the characterization of solutions
of the vector equations.

######  Corollary 2.

If $\mathbf{c}\in\mathcal{C}(\mathbf{A})$, then the general solution of the
vector equation $\mathbf{A}\mathbf{x}=\mathbf{c}$ is

$\displaystyle\mathbf{x}=\mathbf{Y}^{*}\mathbf{c}+(\mathbf{I}-\mathbf{Y^{*}}\mathbf{A})\mathbf{y},$
(49)

where $\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{A})^{+}\mathbf{B}^{*}$ and
$\mathrm{rank}(\mathbf{B}^{*}\mathbf{F})=k$, and $\mathbf{y}$ is arbitrary.

In the following section we show how the meta-factorization procedure works
and how certain well-known matrix factorizations satisfy its equations.

##  3 How meta-factorization describes known factorizations

Meta-factorization explores the structure of known factorizations, thereby
revealing their properties and providing hints for possible modifications.
These insights may become useful for designing numerical algorithms, as they
provide ways to control the properties of numerical algorithms, including
performance, stability, and complexity. The key benefit of meta-factorization
is that it allows for modifying the factors of factorizations while preserving
the ability of the factors’ product to reconstruct the original matrix.

This section illustrates how meta-factorization can be applied to describe and
modify selected factorizations. First, it reconstructs the SVD. Next, it
derives the form of the mixing matrix for the CPQR. Finally, by performing
meta-factorization of the mixing matrix for QR-based decompositions, it
develops a family of UTV-like factorizations.

###  3.1 Singular value decomposition (SVD)

Let us begin by illustrating how meta-factorization procedure works. We apply
it to reconstruct the reduced SVD.

For every rectangular $m\times n$ matrix $\mathbf{A}$, the SVD shows that
$\mathbf{A}$ is isometric to the singular value matrix

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{S}}{\scriptstyle(m\times
n)}=\stackunder{\mathbf{U}^{*}}{\scriptstyle(m\times
m)}&&\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&&\stackunder{\mathbf{V}.}{\scriptstyle(n\times n)}\end{aligned}$ (50)

Matrices $\mathbf{U}$ and $\mathbf{V}$ are complex unitary, and $\mathbf{S}$
is a rectangular diagonal matrix with (real non-negative) singular values on
the diagonal. The fundamental subspaces for $\mathbf{A}$, namely, the column
space, the row space and the corresponding nullspaces, as well as the rank of
$\mathbf{A}$, all emerge as a result of this special factorization.

The reduced singular value decomposition factors a rank-$k$ matrix
$\mathbf{A}$ into the following product:

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&=&\stackunder{\mathbf{U}(:,1{:}k)}{\scriptstyle(m\times
k)}&&\stackunder{\mathbf{S}(1{:}k,1{:}k)}{\scriptstyle(k\times
k)}&&\stackunder{\mathbf{V}(:,1{:}k)^{*}.}{\scriptstyle(k\times
n)}\end{aligned}$ (51)

We can reconstruct the reduced singular value decomposition by following the
steps of meta-factorization. The goal is to build the diagonal mixing matrix
$\mathbf{G}=\mathbf{S}(1{:}k,1{:}k)=\mathrm{diag}(\sigma_{1},\dots,\sigma_{k})$.

The first step is to select bases for $\mathcal{C}(\mathbf{A})$ and
$\mathcal{C}(\mathbf{A}^{*})$. As is well known, the required choice of bases
is

$\displaystyle\mathbf{F}=\mathbf{U}(:,1{:}k)\quad\text{and}\quad\mathbf{H}^{*}=\mathbf{V}(:,1{:}k)^{*},$
(52)

where $\mathbf{U}(:,1{:}k)$ and $\mathbf{V}(:,1{:}k)$ hold orthonormal
eigenvectors corresponding to the nonzero eigenvalues of
$\mathbf{A}\mathbf{A}^{*}$ and $\mathbf{A}^{*}\mathbf{A}$. More generally,
$\mathbf{U}$ and $\mathbf{V}$ are defined by the following (symmetric)
eigenvalue decompositions,

$\displaystyle\begin{aligned}
(\mathbf{A}\mathbf{A}^{*})\mathbf{U}=\mathbf{U}(\mathbf{S}\mathbf{S}^{*})\quad\text{and}\quad(\mathbf{A}^{*}\mathbf{A})\mathbf{V}=\mathbf{V}(\mathbf{S}^{*}\mathbf{S}).\end{aligned}$
(53)

The structure of the SVD decomposition is imposed by the choice of
$\mathbf{F}$ and $\mathbf{H}$. In fact, this particular choice makes
$\mathbf{S}$ diagonal. The columns of
$\mathbf{H}=[\mathbf{v}_{1}\mid\dots\mid\mathbf{v}_{k}]$ are connected with
the columns of $\mathbf{F}=[\mathbf{u}_{1}\mid\dots\mid\mathbf{u}_{k}]$ by the
equations

$\displaystyle\mathbf{A}\mathbf{v}_{i}=\sigma_{i}\mathbf{u}_{i}\quad\text{and}\quad\mathbf{A}^{*}\mathbf{u}_{i}=\sigma_{i}\mathbf{v}_{i}\quad\text{for}\quad
i=1,\dots,k.$ (54)

Thus, by orthonormality of $\mathbf{U}$ (and $\mathbf{V}$), we have

$\displaystyle\mathbf{u}^{*}_{i}\mathbf{A}\mathbf{v}_{i}$
$\displaystyle=\sigma_{i},\ i=1,\dots,k,$ (55)
$\displaystyle\mathbf{u}^{*}_{i}\mathbf{A}\mathbf{v}_{j}$ $\displaystyle=0,\
i\neq j.$ (56)

The second step of meta-factorization should confirm this fundamental fact
and, indeed, it does. Solving the projector equation (5) for $\mathbf{Y}$ and
$\mathbf{X}$, we obtain

$\displaystyle\mathbf{Y}^{*}=\mathbf{U}(:,1{:}k)^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{V}(:,1{:}k).$
(57)

These are precisely the matrices that diagonalize $\mathbf{A}$ and make it
isometric to the diagonal singular value matrix. This is demonstrated in the
final third step of the meta-factorization procedure, which yields

$\displaystyle\begin{aligned}
\mathbf{G}=\mathbf{S}(1{:}k,1{:}k)=\mathbf{U}(:,1{:}k)^{*}\mathbf{A}\mathbf{V}(:,1{:}k)=\mathrm{diag}(\sigma_{1},\dots,\sigma_{k}).\end{aligned}$
(58)

As desired, we have reconstructed the reduced SVD for $\mathbf{A}$, namely,

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{U}(:,1{:}k)\mathbf{U}(:,1{:}k)^{*}\mathbf{A}\mathbf{V}(:,1{:}k)\mathbf{V}(:,1{:}k)^{*}\\\
&=\mathbf{U}(:,1{:}k)\mathbf{S}(1{:}k,1{:}k)\mathbf{V}(:,1{:}k)^{*}.\end{aligned}$
(59)

###  3.2 Column-Pivoted QR decomposition (CPQR)

We now show how meta-factorization reveals the form of the mixing matrix
$\mathbf{G}$ for the CPQR.

The CPQR decomposition factors any rectangular $m\times n$ matrix $\mathbf{A}$
into the product of an orthogonal $m\times m$ matrix $\mathbf{Q}$, an upper
rectangular $m\times n$ matrix $\mathbf{R}$ and a rectangular $n\times n$
permutation matrix $\mathbf{\Pi}$:

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&=&\stackunder{\mathbf{Q}}{\scriptstyle(m\times
m)}&&\stackunder{\mathbf{R}}{\scriptstyle(m\times
n)}&&\stackunder{\mathbf{\Pi}^{*}.}{\scriptstyle(n\times n)}\end{aligned}$
(60)

The CPQR also reveals matrix rank; see also Chan, Gu and Eisenstat [2, 9] for
more details. Given a rank-$k$ matrix $\mathbf{A}$, we obtain

$\displaystyle\begin{aligned}
\mathbf{R}=\begin{bmatrix}\mathbf{R}_{11}&\mathbf{R}_{12}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix},\end{aligned}$ (61)

where $\mathbf{R}_{11}$ is a nonsingular upper triangular $k\times k$ matrix.
The thin QR decomposition is then given by

$\displaystyle\setstackgap{L}{15pt}\begin{aligned}
\stackunder{\mathbf{A}}{\scriptstyle(m\times
n)}&=&\stackunder{\mathbf{Q}(:,1{:}k)}{\scriptstyle(m\times
k)}&&\stackunder{\mathbf{R}(1{:}k,:)}{\scriptstyle(k\times
n)}&&\stackunder{\mathbf{\Pi}^{*}.}{\scriptstyle(n\times n)}\end{aligned}$
(62)

To get additional insight into the structure of the decomposition, we can
perform its meta-factorization. Suppose that $\mathbf{Q}$, $\mathbf{R}$ and
$\mathbf{\Pi}$ are known for $\mathbf{A}$. The goal is to examine the form of
the decomposition’s mixing matrix $\mathbf{G}$.

Define $\mathbf{F}$ and $\mathbf{H}$ by the factors of the thin QR
decomposition,

$\displaystyle\mathbf{F}=\mathbf{Q}(:,1{:}k)\quad\text{and}\quad\mathbf{H}^{*}=\mathbf{R}(1{:}k,:)\mathbf{\Pi}^{*}.$
(63)

Solving the projector equation yields

$\displaystyle\begin{aligned}
\mathbf{G}&=\mathbf{Q}(:,1{:}k)^{*}\mathbf{A}(\mathbf{R}(1{:}k,:)\mathbf{\Pi}^{*})^{+}.\end{aligned}$
(64)

However, by the properties of the thin QR decomposition,

$\displaystyle\begin{aligned}
\mathbf{Q}(:,1{:}k)^{*}\mathbf{A}=\mathbf{R}(1{:}k,:)\mathbf{\Pi}^{*}.\end{aligned}$
(65)

The mixing matrix for the CPQR decomposition is, therefore,

$\displaystyle\mathbf{G}=\mathbf{I}_{k}.$ (66)

Notice that in the light of this result, equation (64) can be viewed as a
special type of diagonalization or similarity-type relation.

By the properties of pseudoinverse, the final step of the meta-factorization
yields the column-pivoted QR decomposition,

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{Q}(:,1{:}k)\mathbf{G}\mathbf{R}(1{:}k,:)\mathbf{\Pi}^{*}\\\
&=\mathbf{Q}(:,1{:}k)\mathbf{R}(1{:}k,:)\mathbf{\Pi}^{*},\end{aligned}$ (67)

as expected.

###  3.3 UTV decompositions

We have already seen how meta-factorization describes selected factorizations,
the SVD and the CPQR. This section shows how it can be applied to design
matrix factorizations. The key step involved is to perform factorization of
the mixing matrix $\mathbf{G}$ encoding the properties of $\mathbf{A}$. We
demonstrate this approach by reconstructing the UTV decomposition and then
introducing its modifications.

As Golub et al. and Stewart presented in [8, 28], matrix $\mathbf{A}$ can be
factored into a product of two unitary matrices, $\mathbf{U}$ and
$\mathbf{V}$, and an upper-triangular or lower-triangular mixing matrix
$\mathbf{T}$. The resulting UTV decomposition222The UTV decomposition with
upper-triangular $\mathbf{T}$ is also known as the URV decomposition, whereas
the one with lower-triangular $\mathbf{T}$ is also known as the ULV
decomposition. is a generalization of the SVD and the CPQR decomposition. The
SVD is the case with diagonal $\mathbf{T}$ and the CPQR is the case with
$\mathbf{V}$ being the permutation matrix.

To reconstruct the UTV decomposition, we only need to consider the row space
projector. Define the row space basis as

$\displaystyle\mathbf{H}^{*}=\mathbf{Q}_{r}(:,1{:}k)^{*}.$ (68)

The projector equation yields

$\displaystyle\mathbf{X}=\mathbf{Q}_{r}(:,1{:}k).$ (69)

As a result, we reach the decomposition

$\displaystyle\mathbf{G}=\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)\quad\text{and}\quad\mathbf{A}=\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)\mathbf{Q}_{r}(:,1{:}k)^{*}.$
(70)

Extending the obtained reconstruction formula for $\mathbf{A}$ to include the
remaining columns of $\mathbf{Q}_{r}$ yields

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)\mathbf{Q}_{r}(:,1{:}k)^{*}\\\
&=\begin{bmatrix}\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)&\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix}\begin{bmatrix}\mathbf{Q}_{r}(:,1{:}k)^{*}\\\
\mathbf{Q}_{r}(:,(k{+}1){:}n)^{*}\end{bmatrix}\\\
&=\begin{bmatrix}\mathbf{G}&\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix}\mathbf{Q}_{r}^{*}.\end{aligned}$
(71)

This is the factorization we wish to transform into the UTV decomposition. The
key step is to perform factorization (meta-factorization) of the mixing matrix
$\mathbf{G}$. In this case, the required factors are provided by the SVD.
Namely, we get

$\displaystyle\begin{aligned}
\mathbf{G}=\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)=\bar{\mathbf{U}}\bar{\mathbf{S}}\bar{\mathbf{V}}^{*}.\end{aligned}$
(72)

Inserting (72) into (71) and factoring out $\bar{\mathbf{U}}$ and
$\bar{\mathbf{V}}^{*}$, we obtain the desired UTV decomposition:

$\displaystyle\begin{aligned}
\mathbf{A}&=\bar{\mathbf{U}}\begin{bmatrix}\bar{\mathbf{S}}&\bar{\mathbf{U}}^{*}\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix}\begin{bmatrix}\bar{\mathbf{V}}^{*}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n-k}\end{bmatrix}\mathbf{Q}_{r}^{*}=\mathbf{U}\mathbf{T}\mathbf{V}^{*},\\\
\mathbf{U}&=\bar{\mathbf{U}},\\\
\mathbf{T}&=\begin{bmatrix}\bar{\mathbf{S}}&\bar{\mathbf{U}}^{*}\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix},\\\
\mathbf{V}&=\mathbf{Q}_{r}\begin{bmatrix}\bar{\mathbf{V}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n-k}\end{bmatrix}.\end{aligned}$ (73)

Notice that $\mathbf{U}$ and $\mathbf{V}$ are indeed unitary (square and
orthonormal), and $\mathbf{T}$ is upper-triangular. For more details on
computing the UTV decomposition see Martinsson et al. [19].

We can now turn to the more general case of factorizations with two sided
projectors. Consider the following CPQR decompositions:

$\displaystyle\mathbf{A}=\mathbf{Q}_{c}\mathbf{R}_{c}\mathbf{\Pi}^{*}_{c}\quad\text{and}\quad\mathbf{A}^{*}=\mathbf{Q}_{r}\mathbf{R}_{r}\mathbf{\Pi}^{*}_{r}.$
(74)

Selecting $\mathbf{Q}_{c}$ to represent the column space and $\mathbf{Q}_{r}$
to represent the row space for $\mathbf{A}$, we have

$\displaystyle\mathbf{F}=\mathbf{Q}_{c}(:,1{:}k)\quad\text{and}\quad\mathbf{H}=\mathbf{Q}_{r}(:,1{:}k),$
(75)

and solving the projector equation, we get

$\displaystyle\mathbf{Y}^{*}=\mathbf{Q}_{c}(:,1{:}k)^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{Q}_{r}(:,1{:}k).$
(76)

From that, we immediately obtain the desired meta-factorization equations,

$\displaystyle\begin{aligned}
\mathbf{G}&=\mathbf{Q}_{c}(:,1{:}k)^{*}\mathbf{A}\mathbf{Q}_{r}(:,1{:}k),\\\
\mathbf{A}&=\mathbf{Q}_{c}(:,1{:}k)\mathbf{Q}_{c}(:,1{:}k)^{*}\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)\mathbf{Q}_{r}(:,1{:}k)^{*}.\end{aligned}$
(77)

Since the SVD-defining condition (54) does not hold in this case, the mixing
matrix $\mathbf{G}$ is not diagonal anymore. However, it does store enough
information to reconstruct $\mathbf{A}$ in the system of coordinates given by
$\mathbf{F}$ and $\mathbf{H}$. This is the observation we can exploit.
Furthermore, we can drop the requirement for both matrices, $\mathbf{U}$ and
$\mathbf{V}$, to be unitary. The result is the design of UTV-like
factorizations.

The key step is factoring the mixing matrix. Consider again performing the SVD
of $\mathbf{G}$:

$\displaystyle\mathbf{G}=\mathbf{Q}_{c}(:,1{:}k)^{*}\mathbf{A}\mathbf{Q}_{r}(:,1{:}k)=\bar{\mathbf{U}}\bar{\mathbf{S}}\bar{\mathbf{V}}^{*}.$
(78)

Introducing that result into (73) defines the following factorization:

$\displaystyle\begin{aligned}
\mathbf{A}&=\bar{\mathbf{U}}\begin{bmatrix}\bar{\mathbf{S}}&\bar{\mathbf{U}}^{*}\mathbf{Q}_{c}(:,1{:}k)^{*}\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix}\begin{bmatrix}\bar{\mathbf{V}}^{*}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n-k}\end{bmatrix}\mathbf{Q}_{r}^{*}=\mathbf{U}\mathbf{T}\mathbf{V}^{*},\\\
\mathbf{U}&=\mathbf{Q}_{c}(:,1{:}k)\bar{\mathbf{U}},\\\
\mathbf{T}&=\begin{bmatrix}\bar{\mathbf{S}}&\bar{\mathbf{U}}^{*}\mathbf{Q}_{c}(:,1{:}k)^{*}\mathbf{A}\mathbf{Q}_{r}(:,(k{+}1){:}n)\end{bmatrix},\\\
\mathbf{V}&=\mathbf{Q}_{r}\begin{bmatrix}\bar{\mathbf{V}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{n-k}\end{bmatrix}.\end{aligned}$ (79)

The following Matlab code illustrates the design:

        F       = Qc(:,1:k);        B       = F;
        H       = Qr(:,1:k);        D       = H;
        [QY,RY] = qr(B’*F,0);       [QX,RX] = qr(H’*D,0);
        Y       = (RY\(QY’*B’))’;   X       = (D/RX)*QX’;
        G       = Y’*A*X;
        [Ubar,Sbar,Vbar] = svd(G);
        U       = Y*Ubar;
        T       = [Sbar,Ubar’*Y’*A*Qr(:,k+1:end)];
        V       = Qr*[Vbar,zeros(k,n-k);zeros(k,n-k)’,eye(n-k)];
        Ar      = U*T*V’;

Notice that $\mathbf{U}$ is not unitary anymore. If $\mathbf{A}$ is a full
column-rank matrix, then $\mathbf{Q}_{c}(:,1{:}k)$ is rectangular and

$\displaystyle\mathbf{U}\mathbf{U}^{*}=\mathbf{Q}_{c}(:,1{:}k)\mathbf{Q}_{c}(:,1{:}k)^{*},$
(80)

which is not an identity matrix in general.

By following the same technique of meta-factorization, we can design another
UTV-like decomposition. Instead of the SVD, performing the CPQR decomposition
of the mixing matrix yields

$\displaystyle\begin{aligned}
\mathbf{G}=\bar{\mathbf{Q}}\bar{\mathbf{R}}\bar{\mathbf{\Pi}}^{*}.\end{aligned}$
(81)

The desired factorization can now be computed as

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{Q}_{c}(:,1{:}k)\bar{\mathbf{Q}}\bar{\mathbf{R}}\bar{\mathbf{\Pi}}^{*}\mathbf{Q}_{r}(:,1{:}k)^{*}=\mathbf{U}\mathbf{T}\mathbf{V}^{*},\\\
\mathbf{U}&=\mathbf{Q}_{c}(:,1{:}k)\bar{\mathbf{Q}},\\\
\mathbf{T}&=\bar{\mathbf{R}},\\\
\mathbf{V}&=\mathbf{Q}_{r}(:,1{:}k)\bar{\mathbf{\Pi}}.\end{aligned}$ (82)

In this case, $\mathbf{U}$ and $\mathbf{V}$ have orthogonal columns and
$\mathbf{T}$ is upper-triangular. The Matlab code illustrates the design:

        F       = Qc(:,1:k);      B       = F;
        H       = Qr(:,1:k);      D       = H;
        [QY,RY] = qr(B’*F,0);     [QX,RX] = qr(H’*D,0);
        Y       = (RY\(QY’*B’))’; X       = (D/RX)*QX’;
        G       = Y’*A*X;
        [Qbar,Rbar,Pbar] = qr(G);
        U       = Qc(:,1:k)*Qbar;
        T       = Rbar;
        V       = Qr(:,1:k)*Pbar;
        Ar      = U*T*V’;

Similarly, to obtain the decomposition with the lower-triangular mixing matrix
$\mathbf{T}$, consider the permuted LU decomposition of the mixing matrix.
This gives rise to the following equality,

$\displaystyle\begin{aligned}
\mathbf{G}=\tilde{\mathbf{\Pi}}^{*}\tilde{\mathbf{L}}\tilde{\mathbf{U}}.\end{aligned}$
(83)

By following the same procedure as before, we obtain

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{Q}_{c}(:,1{:}k)\tilde{\mathbf{\Pi}}^{*}\tilde{\mathbf{L}}\tilde{\mathbf{U}}\mathbf{Q}_{r}(:,1{:}k)^{*}=\mathbf{U}\mathbf{T}\mathbf{V}^{*},\\\
\mathbf{U}&=\mathbf{Q}_{c}(:,1{:}k)\tilde{\mathbf{\Pi}}^{*},\\\
\mathbf{T}&=\bar{\mathbf{L}},\\\
\mathbf{V}&=\mathbf{Q}_{r}(:,1{:}k)\tilde{\mathbf{U}}^{*}.\end{aligned}$ (84)

As before, the Matlab code illustrates the design:

        F       = Qc(:,1:k);      B       = F;
        H       = Qr(:,1:k);      D       = H;
        [QY,RY] = qr(B’*F,0);     [QX,RX] = qr(H’*D,0);
        Y       = (RY\(QY’*B’))’; X       = (D/RX)*QX’;
        G       = Y’*A*X;
        [Lbar,Ubar,Pbar] = lu(G);
        U       = Qc(:,1:k)*Pbar’;
        T       = Lbar;
        V       = Qr(:,1:k)*Ubar’;
        Ar      = U*T*V’;

In that case, $\mathbf{T}$ is lower-triangular by design. However, only
$\mathbf{U}$ has orthogonal columns. That is not the case with $\mathbf{V}$
due to the upper-triangular factor given by $\tilde{\mathbf{U}}$.

This section demonstrated that meta-factorization reconstructs known
factorizations, reveals their internal structures, and allows for introducing
modifications. Let us see if there is more to the application of presented
concept.

##  4 Additional benefits of meta-factorization

This section takes the prospect of meta-factorization to investigate a
relationship of meta-factorization with the theory of generalized matrix
inverses, reconstructs the generalized Nyström’s method for computing
randomized low-rank matrix approximations, and shows the role of the CUR
decomposition.

###  4.1 Pseudoinverse and its formulas

It can be noticed that meta-factorization equations resemble the first and
second Penrose equation characterizing the pseudoinverse of a matrix, which
suggests a relationship of meta-factorization with the theory of generalized
matrix inverses. This section provides some basic observations on this matter.
First, it shows that meta-factorization equations and Penrose equations are
equivalent for square matrices. Second, it shows how meta-factorization can be
used to derive explicit formulas for a pseudoinverse.

The Moore-Penrose pseudoinverse satisfies the two equations

$\displaystyle\begin{aligned}
\mathbf{A}\mathbf{A}^{\\!+}\mathbf{A}=\mathbf{A}\quad\text{and}\quad\mathbf{A}^{\\!+}\mathbf{A}\mathbf{A}^{\\!+}\\!=\mathbf{A}^{\\!+}.\end{aligned}$
(85)

These become meta-factorization equations (under the generalized assumptions
of Theorem 3). Indeed, taking

$\displaystyle\mathbf{F}=\mathbf{A}\quad\text{and}\quad\mathbf{H}^{*}=\mathbf{A}$
(86)

immediately results in the mixing matrix:

$\displaystyle\mathbf{G}=\mathbf{A}^{\\!+}\mathbf{A}\mathbf{A}^{\\!+}\\!=\mathbf{A}^{\\!+},$
(87)

where $\mathbf{Y}^{*}=\mathbf{X}=\mathbf{A}^{\\!+}$. Therefore, we also have

$\displaystyle\mathbf{A}=\mathbf{A}\mathbf{G}\mathbf{A}=\mathbf{A}\mathbf{A}^{\\!+}\mathbf{A}.$
(88)

We can now use meta-factorization to describe the internal structure of
pseudoinverse.

To illustrate the general idea, let us follow James [13], Strang and Drucker
[31], and start with the CR decomposition,

$\displaystyle\mathbf{A}=\mathbf{C}\mathbf{R}.$ (89)

When $\mathbf{A}$ is a rank-$k$ matrix, it can be factored into the product of
the $m\times k$ matrix $\mathbf{C}$, containing the first $k$ independent
columns of $\mathbf{A}$, and the $k\times n$ matrix $\mathbf{R}$, containing
$k$ nonzero rows of the row reduced echelon form of $\mathbf{A}$. Since both
$\mathbf{C}$ and $\mathbf{R}$ have full-rank $k$, the pseudoinverse of
$\mathbf{A}$ is

$\displaystyle\mathbf{A}^{\\!+}=\mathbf{R}^{+}\mathbf{C}^{+}=\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{-1}(\mathbf{C}^{*}\mathbf{C})^{-1}\mathbf{C}^{*}.$
(90)

However, we also have

$\displaystyle(\mathbf{R}\mathbf{R}^{*})^{-1}(\mathbf{C}^{*}\mathbf{C})^{-1}=(\mathbf{C}^{*}\mathbf{C}\mathbf{R}\mathbf{R}^{*})^{-1}=(\mathbf{C}^{*}\mathbf{A}\mathbf{R}^{*})^{-1}.$
(91)

Therefore, we obtain the following formula for the pseudoinverse:

$\displaystyle\mathbf{A}^{\\!+}=\mathbf{R}^{*}(\mathbf{C}^{*}\mathbf{A}\mathbf{R}^{*})^{-1}\mathbf{C}^{*}.$
(92)

We can now derive that formula by following the steps of meta-factorization,
thereby revealing the internal structure of a pseudoinverse. Selecting

$\displaystyle\mathbf{F}=\mathbf{R}^{*}\quad\text{and}\quad\mathbf{H}=\mathbf{C},$
(93)

and solving the projector equation for $\mathbf{Y}$ and $\mathbf{X}$, yields

$\displaystyle\mathbf{Y}^{*}=(\mathbf{R}\mathbf{R}^{*})^{-1}\mathbf{R}\quad\text{and}\quad\mathbf{X}=\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{-1}.$
(94)

This brings us to the decomposition

$\displaystyle\begin{aligned}
\mathbf{G}=(\mathbf{R}\mathbf{R}^{*})^{-1}\mathbf{R}\mathbf{A}^{+}\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{-1}\quad\text{and}\quad\mathbf{A}^{\\!+}=\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{-1}\mathbf{R}\mathbf{A}^{+}\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{-1}\mathbf{C}^{*}.\end{aligned}$
(95)

To show that this is equivalent to (92), it is enough to notice that

$\displaystyle\mathbf{A}^{\\!+}=\mathbf{R}^{+}\mathbf{C}^{+}.$ (96)

Indeed, we have

$\displaystyle\begin{aligned}
\mathbf{A}^{\\!+}&=\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{-1}\mathbf{R}\mathbf{A}^{+}\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{-1}\mathbf{C}^{*}\\\
&=\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{-1}\mathbf{R}\mathbf{R}^{+}\mathbf{C}^{+}\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{-1}\mathbf{C}^{*}\\\
&=\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{-1}(\mathbf{C}^{*}\mathbf{C})^{-1}\mathbf{C}^{*}\\\
&=\mathbf{R}^{*}(\mathbf{C}^{*}\mathbf{A}\mathbf{R}^{*})^{-1}\mathbf{C}^{*},\end{aligned}$
(97)

as required. Thus, the mixing matrix for the pseudoinverse is given by

$\displaystyle\begin{aligned}
\mathbf{G}&=(\mathbf{C}^{*}\mathbf{A}\mathbf{R}^{*})^{-1}.\end{aligned}$ (98)

The following result, presented by Ben-Israel and Greville in [1], shows that
formula (92) is a special case of the general formula derived by MacDuffee.

######  Theorem 7 (MacDuffee [1]).

If a complex $m\times n$ matrix $\mathbf{A}$ of rank $k>0$ has a full-rank
factorization

$\displaystyle\begin{aligned} \mathbf{A}=\mathbf{B}\mathbf{D},\end{aligned}$
(99)

then

$\displaystyle\begin{aligned}
\mathbf{A}^{\\!+}=\mathbf{D}^{*}(\mathbf{B}^{*}\mathbf{A}\mathbf{D}^{*})^{-1}\mathbf{B}^{*}.\end{aligned}$
(100)

Let us derive the formula by following the steps of meta-factorization. By
Theorem 2, orthogonal projections corresponding to

$\displaystyle\mathbf{F}=\mathbf{D}^{*}\quad\text{and}\quad\mathbf{H}^{*}=\mathbf{B}^{*}$
(101)

are given by

$\displaystyle\mathbf{Y}^{*}=(\mathbf{D^{*}})^{+}\quad\text{and}\quad\mathbf{X}=(\mathbf{B}^{*})^{+}.$
(102)

Therefore, the mixing matrix for a pseudoinverse is

$\displaystyle\mathbf{G}=(\mathbf{D}^{*})^{+}\mathbf{A}^{\\!+}(\mathbf{B}^{*})^{+}=(\mathbf{B}^{*}\mathbf{A}\mathbf{D}^{*})^{+}.$
(103)

Since $\mathbf{B}^{*}\mathbf{A}\mathbf{D}^{*}$ is nonsingular as a product of
nonsingular matrices, we have

$\displaystyle\mathbf{A}^{\\!+}$
$\displaystyle=\mathbf{D}^{*}(\mathbf{B}^{*}\mathbf{A}\mathbf{D}^{*})^{-1}\mathbf{B}^{*},$
(104)

giving the desired explicit formula for a pseudoinverse.

###  4.2 Generalized Nyström’s method, the CUR and the outer-product
decomposition

Nakatsukasa [22] showed that many algorithms of randomized linear algebra can
be written in the form of the generalized Nyström’s method. Matrix
$\mathbf{A}$ can be approximated by matrix $\mathbf{A}_{r}$ being a two-sided
projection defined by two sampling (sketch) matrices, an $n\times k$ matrix
$\mathbf{\Omega}_{c}$ and $m\times k$ matrix $\mathbf{\Omega}_{r}$:

$\displaystyle\begin{aligned}
\mathbf{A}_{r}&=\mathbf{A}\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*}\mathbf{A}.\end{aligned}$
(105)

There are two projection matrices involved. The first one,

$\displaystyle\begin{aligned}
\mathbf{P}&=\mathbf{A}\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*},\end{aligned}$
(106)

defines an oblique projection onto the column space of
$\mathbf{A}\mathbf{\Omega}_{c}$. The second,

$\displaystyle\begin{aligned}
\mathbf{R}&=\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*}\mathbf{A},\end{aligned}$
(107)

defines an oblique projection onto the row space of
$\mathbf{\Omega}_{r}^{*}\mathbf{A}$.

The generalized Nyström’s method can be derived from the meta-factorization
procedure. The goal is to design the oblique projectors. By Theorem 2,

$\displaystyle\mathbf{F}=\mathbf{A}\mathbf{\Omega}_{c}\quad$
$\displaystyle\text{and}\quad\mathbf{H}^{*}=\mathbf{\Omega}_{r}^{*}\mathbf{A}$
(108) $\displaystyle\mathbf{B}=\mathbf{\Omega}_{r}\quad$
$\displaystyle\text{and}\quad\mathbf{D}=\mathbf{\Omega}_{c}$ (109)

yield the following solutions to the projector equation:

$\displaystyle\mathbf{Y}^{*}=(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*}\quad\text{and}\quad\mathbf{X}=\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}.$
(110)

As a result, the mixing matrix becomes:

$\displaystyle\begin{aligned}
\mathbf{G}&=(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}=(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+},\end{aligned}$
(111)

which gives rise to the generalized Nyström’s method:

$\displaystyle\begin{aligned}
\mathbf{A}_{r}&=\mathbf{F}\mathbf{G}\mathbf{H}^{*}\\\
&=\mathbf{A}\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{\Omega}_{r}^{*}\mathbf{A}.\end{aligned}$
(112)

As pointed out by Nakatsukasa [22], the preferred way to compute
$\mathbf{A}_{r}$ is to perform QR factorization
$\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c}=\mathbf{Q}\mathbf{R}$
and then take

$\displaystyle\begin{aligned}
\mathbf{A}_{r}&=(\mathbf{A}\mathbf{\Omega}_{c}\mathbf{R}^{-1})(\mathbf{Q}^{*}\mathbf{\Omega}_{r}^{*}\mathbf{A}),\end{aligned}$
(113)

which completes the third step of the meta-factorization procedure. The
approximation error can be substantially reduced with the row-space
oversampling. The Matlab code below illustrates the idea,

        D       = randn(n,k);   B  = randn(m,ceil(2*k)); % row-space oversampling
        F       = A*D;          H  = A’*B;
        [Q,R]   = qr(B’*A*D,0);
        Ar      = (F/R)*(Q’*H’);

The generalized Nyström’s method is closely related to the outer-product
decomposition introduced by Wedderburn [35], Egerváry [7], Householder [12],
Cline and Funderlic [4], and Chu et al. [3]. Suppose that $\mathbf{F}$,
$\mathbf{G}$ and $\mathbf{H}$ are given. Then,

$\displaystyle\begin{aligned}
\mathrm{rank}(\mathbf{A}-\mathbf{F}\mathbf{G}\mathbf{H}^{*})=\mathrm{rank}(\mathbf{A})-\mathrm{rank}(\mathbf{F}\mathbf{G}\mathbf{H}^{*})\end{aligned}$
(114)

if and only if there exist matrices $\mathbf{\Omega}_{c}$ and
$\mathbf{\Omega}_{r}$ such that

$\displaystyle\begin{aligned}
\mathbf{F}=\mathbf{A}\mathbf{\Omega}_{c}\quad\text{and}\quad\mathbf{G}^{-1}=\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c}\quad\text{and}\quad\mathbf{H}=\mathbf{A}^{*}\mathbf{\Omega}_{r}.\end{aligned}$
(115)

Therefore, when the conditions above hold, there exists a finite rank-reducing
process which decomposes $\mathbf{A}$ into a sum of low-rank matrices. The
result is the outer-product decomposition. To be more specific, consider the
Wedderburn rank-one reduction process,

$\displaystyle\begin{aligned} \mathbf{A}_{1}&=\mathbf{A},\\\
\mathbf{A}_{r+1}&=\mathbf{A}_{r}-g_{r}^{-1}\mathbf{A}_{r}\mathbf{u}_{r}\mathbf{v}^{*}_{r}\mathbf{A}_{r},\\\
g_{r}&=\mathbf{v}^{*}_{r}\mathbf{A}_{r}\mathbf{u}_{r}\neq
0\quad\text{for}\quad r=1,\dots,k-1.\end{aligned}$ (116)

The process terminates after $k$ iterations for any rank-$k$ matrix and gives
rise to the factorization of the following form:

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{F}\mathbf{G}\mathbf{H}^{*}=\begin{bmatrix}\mathbf{A}_{1}\mathbf{u}_{1},\dots,\mathbf{A}_{k}\mathbf{u}_{k}\end{bmatrix}\begin{bmatrix}g_{1}&\\\
&\ddots&\\\
&&g_{k}\end{bmatrix}^{-1}\begin{bmatrix}\mathbf{v}^{*}_{1}\mathbf{A}_{1}\\\
\dots\\\
\mathbf{v}^{*}_{k}\mathbf{A}_{k}\end{bmatrix}=\mathbf{A}\mathbf{\Omega}_{c}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{-1}\mathbf{\Omega}_{r}^{*}\mathbf{A}.\end{aligned}$
(117)

We have reached the following conclusion. The generalized Nyström’s method
produces high-quality low-rank approximations of $\mathbf{A}$ by exploiting
the structure of the outer-product factorization and relaxing the rank-
reduction conditions, defined by (115), with random sketch matrices
$\mathbf{\Omega}_{c}$ and $\mathbf{\Omega}_{r}$. In both cases, the factors of
the decompositions satisfy the equations of meta-factorization.

We close the section with one final observation. It regards the relation
between the generalized Nyström’s method and the CUR factorization. Recall
that the CUR factorization provides an approximation:

$\displaystyle\mathbf{A}\approx\mathbf{C}\mathbf{U}\mathbf{R},$ (118)

where $\mathbf{C}=\mathbf{A}(:,J)$ is a subset of representative columns of
$\mathbf{A}$, $\mathbf{R}=\mathbf{A}(I,:)$ is a subset of representative rows
of $\mathbf{A}$, and the mixing matrix is $\mathbf{U}$. The properties of the
approximation depend critically on the rank of $\mathbf{U}$ and the selection
of columns and rows, $J$ and $I$, as demonstrated by Sorensen and Embree [27],
Wang and Zhang [34], Martinsson and Tropp [18], Hamm and Huang [11].

To reconstruct the CUR decomposition from the meta-factorization equations,
consider the following special case of sampling matrices,

$\displaystyle\bar{\mathbf{\Omega}}_{c}=\mathbf{I}_{n}(:,J)\quad\text{and}\quad\bar{\mathbf{\Omega}}_{r}=\mathbf{I}_{m}(:,I).$
(119)

Matrix

$\displaystyle\mathbf{F}=\mathbf{A}\bar{\mathbf{\Omega}}_{c}=\mathbf{C}$ (120)

contains the original columns of $\mathbf{A}$. Similarly, matrix

$\displaystyle\mathbf{H}^{*}=\bar{\mathbf{\Omega}}_{r}^{*}\mathbf{A}=\mathbf{R}$
(121)

contains the original rows of $\mathbf{A}$. For a feasible matrix
$\mathbf{B}$, we then obtain

$\displaystyle\mathbf{Y}^{*}=(\mathbf{B}^{*}\mathbf{A}\bar{\mathbf{\Omega}}_{c})^{+}\mathbf{B}^{*},$
(122)

and for a feasible matrix $\mathbf{D}$, we obtain

$\displaystyle\mathbf{X}=\mathbf{D}(\bar{\mathbf{\Omega}}_{r}^{*}\mathbf{A}\mathbf{D})^{+}.$
(123)

By the rules of meta-factorization, the corresponding mixing matrix is

$\displaystyle\begin{aligned}
\mathbf{G}&=(\mathbf{B}^{*}\mathbf{A}\mathbf{\Omega}_{c})^{+}\mathbf{B}^{*}\mathbf{A}\mathbf{D}(\mathbf{\Omega}_{r}^{*}\mathbf{A}\mathbf{D})^{+}\\\
&=(\mathbf{B}^{*}\mathbf{C})^{+}\mathbf{B}^{*}\mathbf{A}\mathbf{D}(\mathbf{R}\mathbf{D})^{+}.\end{aligned}$
(124)

Selecting

$\displaystyle\mathbf{B}=\mathbf{A}\bar{\mathbf{\Omega}}_{c}=\mathbf{C}\quad$
$\displaystyle\text{and}\quad\mathbf{D}=\mathbf{A}^{*}\bar{\mathbf{\Omega}}_{r}=\mathbf{R}^{*}$
(125)

immediately leads to the following well-known form of the CUR decomposition,

$\displaystyle\begin{aligned}
\mathbf{A}&=\mathbf{C}(\mathbf{C}^{*}\mathbf{C})^{+}\mathbf{C}^{*}\mathbf{A}\mathbf{R}^{*}(\mathbf{R}\mathbf{R}^{*})^{+}\mathbf{R}\\\
&=\mathbf{C}\mathbf{C}^{+}\mathbf{A}\mathbf{R}^{+}\mathbf{R}.\end{aligned}$
(126)

We have reconstructed the CUR decomposition from the meta-factorization
equations.

We can now explain the relation between the CUR factorization and generalized
Nyström’s method. It is defined by the structure of the hidden projectors. The
projectors defined by (125) give rise to the CUR decomposition, whereas the
projectors defined by (109) give rise to the generalized Nyström’s method.

The following Matlab code illustrates the concepts for a random (naive)
selection of columns and rows:

        In      = eye(n);       Im      = eye(m);
        I       = randperm(m,k);J       = randperm(n,k);
        Q       = In(:,J);      P       = Im(:,I);
        F       = A*Q;          H       = A’*P;
        B       = P;            D       = Q;              % Nystrom-like method, or
        B       = F;            D       = H;              % CUR (with naive sampling)
        [QY,RY] = qr(B’*F,0);   [QX,RX] = qr(H’*D,0);     % If rank(B’*F) = rank(H’*D) = k,
        Y       = (RY\(QY’*B’))’;   X   = (D/RX)*QX’;     %  calculate Y and X.
        Y       = (pinv(B’*F)*B’)’; X   = D*pinv(H’*D);   % Otherwise, consider using pinv()
        G       = Y’*A*X;
        Ar      = F*G*H’;

The examples of this section provide valuable insights into the world of
matrix factorizations and reveal relations between well-known matrix
structures and randomized algorithms.

##  5 Summary

The concept of meta-factorization develops a perspective in which matrix
factorizations become solutions of linear matrix equations. This perspective
leads to exciting insights when applied in the study of matrix algorithms. It
allows for reconstructing known factorizations, investigating their internal
structures, and introducing modifications on demand. Furthermore, once the
prospect of meta-factorization is considered, it is possible to notice
analogies between theories previously hidden from view. That is how we have
gained insights investigating the generalized matrix inverses and randomized
linear algebra algorithms. Also, there is much more to be studied about that
idea and further twists are yet to be discovered.

## Acknowledgment

This paper was written in the tumultuous times of the COVID-19 pandemic, where
many plans had to be changed. Nevertheless, I would like to thank Prof.
Gilbert Strang for his enthusiasm and helpful remarks during our remote
discussions that covered my postponed stay at MIT. I also thank Prof. Michael
Saunders for his support and review of the manuscript, Prof. Yuji Nakatsukasa
for helpful discussions, and Prof. Alex Townsend for important references.
Finally, this paper would not have been written without the support of my
wife, Dr. Inez Okulska, who helped me turn a jungle of thoughts into a paper
that hopefully reads well despite its terrifying amount of equations.

##  Appendix A Side note on generalizations

The projector equation plays one of the leading roles in the developed theory
of meta-factorization. Therefore, it is tempting to ask if it can be
generalized and what would be the outcomes of such a generalization? One
potential answer is included here, possibly for future investigation.

Consider a process of recursive meta-factorization defined by a composition of
meta-factorization equations. For the fixed column space and row space bases,
$\mathbf{F}$ and $\mathbf{H}^{*}$, the reconstruction formula for $\mathbf{A}$
becomes

$\displaystyle\begin{aligned}
\mathbf{A}&=(\mathbf{F}\mathbf{Y}^{*})(\mathbf{F}\mathbf{Y}^{*})\cdots(\mathbf{F}\mathbf{Y}^{*})\mathbf{A}(\mathbf{X}\mathbf{H}^{*})\cdots(\mathbf{X}\mathbf{H}^{*})(\mathbf{X}\mathbf{H}^{*})\\\
&=(\mathbf{F}\mathbf{Y}^{*})^{N}\mathbf{A}(\mathbf{X}\mathbf{H}^{*})^{N}\\\
&=\mathbf{F}(\mathbf{Y}^{*}\mathbf{F})^{N-1}\mathbf{Y}^{*}\mathbf{A}\mathbf{X}(\mathbf{H}^{*}\mathbf{X})^{N-1}\mathbf{H}^{*}.\end{aligned}$
(127)

Suppose that we generalize the projector equation, so that we have:

$\displaystyle\begin{aligned}
\mathbf{Y}^{*}\mathbf{F}=\mathbf{Z}_{c}\quad\text{and}\quad\mathbf{H}^{*}\mathbf{X}=\mathbf{Z}_{r}\quad\text{and}\quad\mathbf{Z}_{c}^{N}=\mathbf{Z}_{r}^{N}=\mathbf{I}_{k}.\end{aligned}$
(128)

Observe now that, given rank-$k$ matrices $\mathbf{F}$ and $\mathbf{H}^{*}$,
the generalized projector equations admit the following solutions:

$\displaystyle\begin{aligned}
\mathbf{Y}^{*}=\mathbf{Z}_{c}\mathbf{F}^{+}\quad\text{and}\quad\mathbf{X}=(\mathbf{H}^{*})^{+}\mathbf{Z}_{r}.\end{aligned}$
(129)

We can now see that as a result we get:

$\displaystyle\begin{aligned}
(\mathbf{F}\mathbf{Y}^{*})^{N}=\mathbf{F}\mathbf{Z}_{c}^{N-1}\mathbf{Y}^{*}=\mathbf{F}\mathbf{Z}_{c}^{N}\mathbf{F}^{+}=\mathbf{F}\mathbf{F}^{+}\end{aligned}$
(130)

and:

$\displaystyle\begin{aligned}
(\mathbf{X}\mathbf{H}^{*})^{N}=\mathbf{X}\mathbf{Z}_{r}^{N-1}\mathbf{H}^{*}=(\mathbf{H}^{*})^{+}\mathbf{Z}_{r}^{N}\mathbf{H}^{*}=(\mathbf{H}^{*})^{+}\mathbf{H}^{*}.\end{aligned}$
(131)

In other words, for $p=0,1,2,\dots$ we obtain:

$\displaystyle\begin{aligned}
\mathbf{A}=(\mathbf{F}\mathbf{Y}^{*})^{Np}\mathbf{A}(\mathbf{X}\mathbf{H}^{*})^{Np}=\mathbf{F}\mathbf{F}^{+}\mathbf{A}(\mathbf{H}^{*})^{+}\mathbf{H}^{*}.\end{aligned}$
(132)

Therefore, any matrix factorization:

$\displaystyle\begin{aligned}
\mathbf{A}=\mathbf{F}\mathbf{G}\mathbf{H}^{*}\end{aligned}$ (133)

becomes periodic (with period $N$), when:

$\displaystyle\begin{aligned}
\mathbf{G}=\mathbf{Z}_{c}\mathbf{F}^{+}\mathbf{A}(\mathbf{H}^{*})^{+}\mathbf{Z}_{r}\quad\text{and}\quad\mathbf{Z}_{c}^{N}=\mathbf{Z}_{r}^{N}=\mathbf{I}_{k}.\end{aligned}$
(134)

This is one of the potential spots to investigate the internal structure of
new factorizations333For more details and numerical demonstration see:
https://github.com/mkarpowi/mft. Meta-factorization provides hints on the
structure and offers the potential of creating new ones.

## References

  * [1] Adi Ben-Israel and Thomas N.E. Greville.  Generalized inverses: theory and applications.  Springer Science & Business Media, 2003.
  * [2] Tony F. Chan.  Rank-revealing QR factorizations.  Linear algebra and its applications, 88:67–82, 1987.
  * [3] Moody T Chu, Robert E Funderlic, and Gene H Golub.  A rank–one reduction formula and its applications to matrix factorizations.  SIAM review, 37(4):512–530, 1995.
  * [4] Randall E Cline and Robert E Funderlic.  The rank of a difference of matrices and associated generalized inverses.  Linear Algebra and its Applications, 24:185–215, 1979.
  * [5] James W. Demmel.  Applied Numerical Linear Algebra.  SIAM, 1997.
  * [6] Alan Edelman and Sungwoo Jeong.  Fifty three matrix factorizations: A systematic approach.  arXiv:2104.08669, 2021.
  * [7] Eugen Egerváry.  On rank-diminishing operations and their applications to the solution of linear equations.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 11(5):376–386, 1960.
  * [8] Gene H Golub and Charles F Van Loan.  Matrix Computations.  JHU Press, 2013.
  * [9] Ming Gu and Stanley C Eisenstat.  Efficient algorithms for computing a strong rank-revealing QR factorization.  SIAM Journal on Scientific Computing, 17(4):848–869, 1996.
  * [10] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM Review, 53(2):217–288, 2011.
  * [11] Keaton Hamm and Longxiu Huang.  Perspectives on CUR decompositions.  Applied and Computational Harmonic Analysis, 48(3):1088–1099, 2020\.
  * [12] Alston Scott Householder.  The Theory of Matrices in Numerical Analysis.  Blaisdell, 1964.
  * [13] M James.  The generalised inverse.  The Mathematical Gazette, 62(420):109–114, 1978.
  * [14] Ravindran Kannan and Santosh Vempala.  Randomized algorithms in numerical linear algebra.  Acta Numerica, 26:95, 2017.
  * [15] Carl E. Langenhop.  On generalized inverses of matrices.  SIAM Journal on Applied Mathematics, 15(5):1239–1246, 1967.
  * [16] Michael W Mahoney.  Lecture notes on randomized linear algebra.  arXiv:1608.04481, 2016.
  * [17] Michael W Mahoney and Petros Drineas.  CUR matrix decompositions for improved data analysis.  Proceedings of the National Academy of Sciences, 106(3):697–702, 2009.
  * [18] Per-Gunnar Martinsson.  Randomized methods for matrix computations.  The Mathematics of Data, 25:187–231, 2019.
  * [19] Per-Gunnar Martinsson, Gregorio Quintana-Orti, and Nathan Heavner.  randUTV: A blocked randomized algorithm for computing a rank-revealing UTV factorization.  ACM Transactions on Mathematical Software, 45(1):1–26, 2019.
  * [20] Per-Gunnar Martinsson and Joel Tropp.  Randomized numerical linear algebra: Foundations & algorithms.  arXiv:2002.01387, 2020.
  * [21] George Matsaglia and George PH Styan.  Equalities and inequalities for ranks of matrices.  Linear and multilinear Algebra, 2(3):269–292, 1974.
  * [22] Yuji Nakatsukasa.  Fast and stable randomized low-rank matrix approximation.  arXiv:2009.11392, 2020.
  * [23] Yuji Nakatsukasa and Joel A. Tropp.  Fast & accurate randomized algorithms for linear systems and eigenvalue problems.  arXiv:2111.00113, 2021.
  * [24] Roger Penrose.  A generalized inverse for matrices.  In Mathematical proceedings of the Cambridge Philosophical Society, volume 51, pages 406–413. Cambridge University Press, 1955.
  * [25] Roger Penrose.  On best approximate solutions of linear matrix equations.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 52, pages 17–19. Cambridge University Press, 1956.
  * [26] Gil Shabat, Yaniv Shmueli, Yariv Aizenbud, and Amir Averbuch.  Randomized LU decomposition.  Applied and Computational Harmonic Analysis, 44(2):246–272, 2018\.
  * [27] Danny C. Sorensen and Mark Embree.  A DEIM induced CUR factorization.  SIAM Journal on Scientific Computing, 38(3):A1454–A1482, 2016.
  * [28] G. W. Stewart.  Matrix Algorithms: Basic Decompositions.  SIAM, 1998.
  * [29] Gilbert Strang.  Computational Science and Engineering.  Wellesley-Cambridge Press, Wellesley, MA, 2007.
  * [30] Gilbert Strang.  Linear Algebra and Learning from Data.  Wellesley-Cambridge Press Cambridge, 2019.
  * [31] Gilbert Strang and Daniel Drucker.  Three matrix factorizations from the steps of elimination.  In preparation, 2021.
  * [32] Alex Townsend and Lloyd N. Trefethen.  Continuous analogues of matrix factorizations.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2173), 2015.
  * [33] Lloyd N. Trefethen and David Bau III.  Numerical Linear Algebra.  SIAM, 1997.
  * [34] Shusen Wang and Zhihua Zhang.  Improving CUR matrix decomposition and the Nyström approximation via adaptive sampling.  The Journal of Machine Learning Research, 14(1):2729–2769, 2013\.
  * [35] Joseph Henry Maclagan Wedderburn.  Lectures on matrices, volume 17.  American Mathematical Soc., 1934.