\begin{document}

\title{\bf  Random Multipliers Numerically Stabilize Gaussian and Block Gaussian Elimination: 
Proofs and an Extension to Low-rank Approximation
\thanks {Some results of this paper have been presented at the 
 ACM-SIGSAM International 
Symposium on Symbolic and Algebraic Computation (ISSAC '2011), in San Jose, CA, 
June 8--11, 2011,
the
3nd International Conference on Matrix Methods in Mathematics and 
Applications (MMMA 2011), in
Moscow, Russia, June 22--25, 2011, 
the 7th International Congress on Industrial and Applied Mathematics 
(ICIAM 2011), in Vancouver, British Columbia, Canada, July 18--22, 
2011,
the SIAM International Conference on Linear Algebra,
in Valencia, Spain, June 18-22, 2012,  
the Conference on Structured Linear and Multilinear Algebra Problems (SLA2012),
in  Leuven, Belgium, September 10-14, 2012, and the 18th Conference of the 
International Linear Algebra Society (ILAS 2013), in Providence, RI, USA, June 3--7, 2013}}

\author{Victor Y. Pan$^{[1, 2],[a]}$, Guoliang Qian$^{[2],[b]}$ and Xiaodong Yan $^{[2],[c]}$


\and\\
$^{[1]}$ Department of Mathematics and Computer Science \\
Lehman College of the City University of New York \\
Bronx, NY 10468 USA \\
$^{[2]}$ Ph.D. Programs in Mathematics  and Computer Science \\
The Graduate Center of the City University of New York \\
New York, NY 10036 USA \\
$^{[a]}$ victor.pan@lehman.cuny.edu \\
http://comet.lehman.cuny.edu/vpan/  \\
$^{[b]}$ gqian@gc.cuny.edu \\
$^{[c]}$ xyannyc@yahoo.com \\
} 


\maketitle


\begin{abstract}

We study two applications of standard Gaussian random multipliers. At 
first we prove that with a probability close to 1 such a multiplier is 
expected to numerically stabilize Gaussian elimination with no pivoting
as well as block Gaussian elimination. Then, by extending our analysis, 
we prove that such a multiplier is also expected to support low-rank 
approximation of a matrix without customary oversampling. Our test 
results are in good accordance with this formal study. The results 
remain similar when we replace Gaussian multipliers with random circulant 
or Toeplitz multipliers, which involve fewer random parameters and enable 
faster multiplication. We formally support the observed efficiency of random 
structured multipliers  applied to approximation, but not to elimination. 
Moreover, we prove that with a probability close to 1 Gaussian random 
circulant multipliers do not fix numerical instability of the elimination 
algorithms for a specific narrow class of well-conditioned inputs. We know 
of no such hard input classes for various alternative choices of random 
structured multipliers, but for none of such multipliers
we have a formal proof of its efficiency  
 for numerical Gaussian elimination.
\end{abstract}

\paragraph{\bf 2000 Math. Subject Classification:}
 15A52, 15A12, 15A06, 65F22, 65F05

\paragraph{\bf Key Words:}
Random matrices,
Random multipliers,
Gaussian elimination,
Pivoting,
Block Gaussian elimination,
Low-rank approximation,
SRFT matrices,
Random circulant matrices

 


\section{Introduction}\label{sintro}


\subsection{Overview}\label{sover}


We call a standard Gaussian random matrix
 just {\em Gaussian}
and  
apply Gaussian multipliers to

\begin{itemize}

\item 
 numerically stabilize
 Gaussian 
and block Gaussian elimination
 using no pivoting, orthogonalization
or symmetrization
\item 
approximate the leading singular spaces
of an ill-con\-di\-tioned matrix, 
 associated with 
its largest singular values 
 and 
\item 
approximate this matrix by a low-rank matrix.
\end{itemize}


Ample empirical evidence shows efficiency of all these applications
 \cite{HMT11}, \cite{M11}, \cite{PQZ13}, 
but formal proofs are only known
in the case of approximation  
under the assumption of small oversampling,
that is, a minor increase of the multiplier size  \cite{HMT11}.
We provide 
formal support for randomized stabilization of Gaussian and block Gaussian 
elimination  and for
randomized approximation without oversampling. 

Our test results are in good accordance with our analysis
and also show similar power of  random circulant
multipliers or their leading (that is, northeastern) Toeplitz blocks 
for the   
elimination and approximation, respectively.  
Simple extension of the known 
estimates enables formal support for this empirical observation
in the case of 
low-rank approximation, but not
in the case of elimination.
Moreover, we prove that with a high probability 
Gaussian circulant multipliers cannot fix numerical instability
 of the elimination algorithms 
for a specific narrow class of inputs.
The issue remains open for the elimination preprocessed 
with other random structured multipliers.

In the next subsections we specify further the computational tasks and our results.


\subsection{Numerical Gaussian elimination with no pivoting and 
block Gaussian elimination}\label{sgenp}


 
Gaussian elimination, applied
numerically,
with rounding errors, can 
  fail even in the case of 
a nonsingular well-con\-di\-tioned input matrix
unless
this matrix is also
positive definite,
 diagonally 
dominant,  or 
totally positive.
For example (see part (i) of Theorem \ref{thcgenp}), 
numerical Gaussian elimination fails when it is 
applied to
the  matrices of discrete Fourier transform of large sizes,
even though they  are unitary up to scaling.
In practice the
user avoids the problems by applying 
 Gaussian elimination 
with partial pivoting,
that is, appropriate row interchange,
which has some limited 
formal and ample empirical support. 
Our alternative is  
Gaussian elimination with no pivoting.
Hereafter we use the acronyms  {\em GEPP} and 
{\em  GENP}.

In Section \ref{ssgnp} we prove that 
GENP is safe numerically
with a probability close to 1
 when the input matrix
 is  
nonsingular, well-conditioned and preprocessed
with a Gaussian 

 multiplier. 

At the first glance our preprocessing may seem to be more costly than
partial pivoting. The latter only involves order of $n^2$
comparisons for an $n\times n$ input matrix $A$, versus 
$2n^3-n^2$ arithmetic operations for its multiplication 
by a Gaussian multiplier.
Generation of $n^2$ independent Gaussian 
entries of the  $n\times n$ multiplier
is  an additional charge, even though 
 one can discount partly it because this  
preprocessing stage is independent of the input matrix.
More careful comparison, however, shows that
partial pivoting takes quite a heavy toll.


It 
interrupts the stream of arithmetic operations 
with foreign operations of comparison,
involves book-keeping, compromises data locality, and
increases communication overhead and data dependence. 

 
Choosing 
between GEPP and GENP with randomized preprocessing 


the user may also consider various other factors,  
 some of them dynamic in time.
For example, here is a relevant citation from  \cite{BCD14}:
``The traditional metric for the efficiency of a numerical algorithm has been the number of arithmetic operations it performs. Technological trends have long been reducing the time to perform an arithmetic operation, so it is no longer the bottleneck in many algorithms; rather, communication, or moving data, is the bottleneck".

Our  test results with Gaussian input  matrices
are in good accordance with our formal estimates
 (see Figures \ref{fig:fig1}--\ref{fig:fig5} 
 and Tables \ref{tab61}--\ref{tab65a}). 


In our tests  
the output accuracy of GENP with preprocessing
was a little lower than in the case of the customary GEPP,
as can be expected,
but 
 a single step of 
iterative refinement,
performed at a dominated 
computational cost, has always fixed this discrepancy
(see Figures \ref{fig:fig2} and \ref{fig:fig3} and  Table \ref{tab63}). 
The tests show similar results when we applied Gaussian circulant
rather than general Gaussian
multipliers, and in this case the cost of our preprocessing 
decreases dramatically
(see more on that in Section \ref{srsm}).

Finally, 
all our study of GENP,  including formal 
probabilistic support of its numerical performance,
is immediately extended to {\em block Gaussian elimination}
(see Section \ref{sbgegenpn}),
whereas  pivoting cannot amend this valuable algorithm,
and so without our preprocessing, it is numerically unsafe to use it,
unless an  input matrix is nonsingular, well-con\-di\-tioned, and
positive definite,
 diagonally 
dominant,  or 
totally positive.
 


\subsection{Low-rank approximation of a matrix}\label{srdpr}


Random multipliers are known to be highly
efficient  
 for low-rank approximations
of an $m\times n$ matrix $A$ having a small numerical rank $r$.
As the basic step, one computes the product $AH$ where $H$ is 
a random $n\times l$ multiplier, for  $l=r+p$ and 
a positive {\em oversampling integer} $p$.
The resulting randomized algorithms  have been studied
extensively, both formally and experimentally
 \cite{HMT11},  \cite{M11}.
They
are
numerically stable, run at a low 
computational cost, 
allow low-cost improvement of 
the output accuracy by means of the Power Method,
and have 
 important applications to
matrix computations, data mining,
statistics, PDEs and integral equations.

By extending our analysis  of preprocessed GENP, 
we prove that even for an $n\times r$ Gaussian
multiplier, that is, {\em without oversampling},
the algorithms output 
rank-$r$  approximations of the input matrix
with a probability close to 1.
Then again our test results 
are in good accordance with our formal estimates
 (see Figures \ref{fig:fig6} 
and \ref{fig:fig7} and Tables \ref{tab67}--\ref{tab614}).

The decrease of $p$ to 0 should be theoretically interesting,
although it has only  minor practical promise, 
limited to the cases where   
the numerical rank $r$ is small and  
the user knows it. Indeed, according to \cite[Section 4.2]{HMT11},
``it is adequate to choose  $\dots~~p=5$ or $p=10$".

\subsection{Computations with random structured multipliers}\label{srsm}

The SRFT $n\times l$  multipliers $H$ 
(SRFT is the acronym for Subsample Random Fourier Transform)
involve only $n$ random parameters
versus $nl$ parameters of Gaussian multipliers and accelerate 
the computation of the product $AH$ by a factor of $l/\log (l)$
 versus $n\times l$ Gaussian multipliers $H$.
It has been proved that they are expected to support rank-$r$ approximation
 assuming
oversampling integers $p=l-r$ of order $r\log (r)$,
and empirically this has been observed for reasonably small constants $p$, usually 
being not more than 20
 \cite[Section 11]{HMT11}, \cite{M11}.
We readily extend the latter results
 to the case where $n\times l$ products of random  $n\times n$
circulant and random $n\times l$ permutation matrices
are used as multipliers instead of SRFT matrices (see Remark \ref{remsrtft}),
and we observe such empirical behavior also in the cases when we 
preprocess GENP by applying random circulant  multipliers
(see  Figures  \ref{fig:fig2},  \ref{fig:fig3}, \ref{fig:fig6} 
and \ref{fig:fig7} and Tables \ref{tab64}, \ref{tab65}, \ref{tab611}, and \ref{tab614}).
In the latter case we only need $n$
random parameters versus $n^2$ for a Gaussian multiplier
and accelerate Gaussian preprocessing  for GENP by a factor of $n/\log (n)$.
  
This acceleration factor grows to $n^2/\log (n)$ 
if an input matrix is Toeplitz or Toeplitz-like
and if we can apply the
 MBA  celebrated 
 algorithm, which is just recursive block
 Gaussian elimination adjusted to a Toeplitz-like input \cite[Chapter 5]{p01}
and which is superfast, that is, runs in nearly linear arithmetic time.
Numerical stability problems, however, are well known for this algorithm
\cite{B85}, and fixing them was the central subject of the highly recognized 
papers \cite{GKO95} and \cite{G98}. Pivoting could not be applied here
because it destroys Toeplitz structure, thus increasing the solution cost to cubic.
So the authors first reduced the task to the  case of Cauchy-like inputs by
specializing the techniques of the transformation of matrix structures from \cite{P90} 
(also see \cite{P15})
and then applied a fast Cauchy-like variant of GEPP
using quadratic arithmetic time. 
 
Stabilization of the MBA superfast algorithm by means  of random
circulant multipliers would mean randomized acceleration 
by order of magnitude versus   \cite{GKO95} and \cite{G98}
because circulant multipliers preserve Toeplitz-like structure
of an input matrix and thus keep the MBA algorithm superfast.
Empirically, Gaussian circulant multipliers do stabilize
GENP numerically,
but no formal support is known for  this observed behavior.
Moreover we
even prove that 
the application  of GENP 
 fails
 for a specific narrow class of  matrices, 
and with probability near 1
 their preprocessing with Gaussian circulant multipliers 
does not help 
 (see Remark \ref{remsrtft}). 
We know of  no such hard inputs
 for some variations 
of using Gaussian circulant multipliers for GENP and MBA, such as
using products of random circulant and skew-circulant multipliers
and using random circulant pre-multipliers and post-multipliers simultaneously.
Both variations preserve Toeplitz-like input structure 
and allow superfast performance of the MBA algorithm,
but we have no formal support for the efficiency of such multipliers.

\subsection{Related works}\label{srel}

Preconditioning of linear systems of equations  
is a classical subject \cite{A94}, \cite{B02}, \cite{G97}.
For early work on randomized multiplicative preprocessing
as a means of countering degeneracy of matrices
see 
Section 2.13 ``Regularization of a Matrix via 
Preconditioning with Randomization" in \cite{BP94}
 and the bibliography therein. On the most recent advance in this direction see \cite{PZa}.
On the early specialization of such techniques to Gaussian elimination
see \cite{PP95}.
Randomized multiplicative preconditioning
for numerical stabilization of GENP 
was proposed in \cite[Section 12.2]{PGMQ} 
and \cite{PQZ13}, although only
weaker theorems on the
 formal support of this approach
were stated and their proofs were omitted.
The paper \cite{BBD12}
and the bibliography therein cover the
heuristic application of PRBMs
(that is, Partial Random Butterfly Multipliers), providing 
some empirical support 
 for GENP with  preprocessing.
On low-rank approximation
we refer the reader to the surveys \cite{HMT11} and \cite{M11},
which were the springboard for our study in Section \ref{sapsr}.
We cite these and other related works throughout 
the paper and refer the reader to \cite[Section 11]{PQZa}
for further bibliography. 
The estimates of our Corollary \ref{cor1} 
are close to the ones of  \cite[Theorem 3.8]{PQ10},
which
were 
the basis for our algorithms in \cite{PQ10}, 
\cite{PQ12}, and \cite{PQZC}.  
Unlike the latter papers, however, we 
state these basic estimates in a simpler
form, refine them by following 
 \cite{CD05} rather than \cite{SST06},
and include their detailed proofs.
On the related subject of estimating the norms and condition
numbers of Gaussian matrices and random structured matrices 
see   \cite{D88},  \cite{E88},  \cite{ES05},
\cite{CD05}, \cite{SST06}, \cite{HMT11}, \cite{T11}, and \cite{PSZa}. 
For a natural extension
of our present work,
 one can combine
randomized matrix multiplication with
randomized augmentation and 
 additive preprocessing
of   \cite{PQ10}, \cite{PQ12}, and
 \cite{PQZC}. 


\subsection{Organization of the paper}\label{sorg}


In the next 
section
we recall some definitions and basic results.
In Section \ref{sbgegenpn} we show that 
GENP and block Gaussian elimination
are numerically safe for a matrix whose 
all leading blocks are nonsingular 
and well conditioned. 
In Section \ref{smrc} we  estimate  
the impact of  
 preprocessing
with general nonrandom multipliers on 
these properties of the leading blocks.
In Sections \ref{sapsr1}
and
\ref{srnd} we extend our analysis 
from  Section \ref{smrc} in order to cover the impact of 
Gaussian and random structured multipliers,
respectively.


In Section \ref{sapsr}
 we recall  an algorithm from  \cite{HMT11} for
low-rank approximation and prove that this 
randomized algorithm is expected to
work even with no oversampling.
 In Section \ref{sexp} we cover numerical tests

(the contribution of the last two authors).
Section \ref{sconcl} contains a brief summary.


In  Appendix \ref{sssm} we recall the known
probabilistic estimates for the error norms of 
randomized low-rank approximations.

In  Appendix \ref{srsnrm}
we estimate the probability that 
 a 
 random matrix 
has full rank
under the uniform probability distribution.
In Appendix \ref{sosvdi0}
we  estimate the perturbation errors of matrix inversion.
In Appendix \ref{stab} we display tables with our test results,
which are more detailed than the data given by the plots in Section  \ref{sexp}.  
Some readers may be only interested in 
the part of our paper on GENP.
They can skip Sections 
 \ref{sapsr}
 and \ref{stails}.


\section{Some definitions}\label{sdef}


Except for using unitary circulant matrices in Sections  \ref{srnd} and
\ref{stails}, we
assume computations in the field $\mathbb R$ of real numbers,
but the extension to the case of the complex field $\mathbb C$
is quite straightforward.

Hereafter ``flop" stands for ``arithmetic operation",
`` i.i.d." stands for ``independent identically distributed",
and ``Gaussian matrix" stands for
``standard Gaussian random matrix" (cf. Definition \ref{defrndm}). 
 The concepts ``large", ``small", ``near", ``closely approximate", 
``ill-con\-di\-tioned" and ``well-con\-di\-tioned" are 
quantified in the context. By saying ``expect" and ``likely" we
mean ``with probability $1$ or close to $1$".
(We only use the concept of the expected value in Theorem \ref{thszcd}, 
Corollary \ref{cocd}, and Appendix \ref{sssm}.) 


Next we recall and extend some customary definitions of matrix computations
\cite{GL13}, \cite{S98}.


$\mathbb  R^{m\times n}$ is the class of real $m\times n$ matrices $A=(a_{i,j})_{i,j}^{m,n}$.


$\diag (B_1,\dots,B_k)=\diag(B_j)_{j=1}^k$ is a $k\times k$ block diagonal matrix 
with the diagonal blocks $B_1,\dots,B_k$. $(B_1)_{1}^k$, 
$(B_1~|~\dots~|~B_k)$, and $(B_1,\dots,B_k)$ 
denote a $1\times k$ block matrix with the blocks $B_1,\dots,B_k$.
In both cases the blocks $B_j$ can be rectangular.
  
$I_n$
is the $n\times n$ identity
matrix. 
$O_{k,l}$ is the $k\times l$ matrix filled with zeros. 
We write $I$ and $O$ if the  matrix size  
 is defined by context. 


$A^T$ is the transpose  of a 
matrix $A$.

 $A_{k,l}$  denotes its leading, 
that is, northwestern  $k\times l$ block 
submatrix, and 

we also write  $A^{(k)}=A_{k,k}$.


$||A||=||A||_2$
is the spectral norm of a matrix $A$.
$||A||_F$
is its Frobenius norm.

A real matrix $Q$ is  
{\em orthogonal} if $Q^TQ=I$ 
 or $QQ^T=I$. 
$(Q,R)=(Q(A),R(A))$ for an $m\times n$ matrix $A$ of rank $n$
denotes a unique pair of orthogonal $m\times n$ and
 upper triangular $n\times n$ matrices such that $A=QR$
and all diagonal entries of the matrix $R$
are positive  \cite[Theorem 5.2.3]{GL13}. 

 $A^+$ denotes the Moore--Penrose 
pseudo-inverse of an $m\times n$ matrix $A$, and
\begin{equation}\label{eqsvd}
A=S_A\Sigma_AT_A^T
\end{equation}
denotes its SVD where $S_A^TS_A=S_AS_A^T=I_m$,
$T_A^TT_A=T_AT_A^T=I_n$,
$\Sigma_A=\diag(\sigma_j(A))_{j}$,  and

$\sigma_j=\sigma_j(A)$

is the $j$th largest singular value of $A$.
If a matrix $A$ has full column rank $\rho$, then


\begin{equation}\label{eqnrm+}
||A^+||=1/\sigma_{\rho}(A).
\end{equation}

  $A^{+T}$ stands for $(A^+)^T=(A^T)^+$, $A_s^T$  for $(A_s)^{T}$,
and $A_s^+$ for $(A_s)^{+}$ where $s$ can denote a scalar, a matrix, or
a pair of such objects, e.g., $A_{k,l}^T$ stands for $(A_{k,l})^T$. 


$\kappa (A)=\frac{\sigma_1(A)}{\sigma_{\rho}(A)}=||A||~||A^+||$ is the condition 
number of an $m\times n$ matrix $A$ of a rank $\rho$. Such matrix is 
{\em ill-con\-di\-tioned}
if the ratio $\sigma_1(A)/\sigma_{\rho}(A)$
is large. If the ratio is reasonably bounded,
then the matrix is {\em well-con\-di\-tioned}.
An $m\times n$ matrix $A$ has a {\em numerical rank} 
$r=\nrank(A)\le \rho=\rank (A)$ 
if the ratios $\sigma_{j}(A)/||A||$
are small for $j>r$ but not for $j\le r$. 


The following concepts 
cover all rectangular matrices, but we need 
them just in the case of square matrices, whose sets of leading blocks 
include the   matrices themselves.
A matrix is {\em strongly nonsingular} if all its leading blocks
 are nonsingular. 
Such a matrix is 
 {\em strongly well-conditioned} if all its leading blocks
 are well-conditioned. 

We recall further relevant definitions and basic results  of
 matrix computations
in the beginning of Section \ref{sapsr} and in the Appendix.


\section{Block Gaussian elimination and GENP
}\label{sbgegenpn}


For a 
nonsingular
$2\times 2$ block matrix $A=\begin{pmatrix}
B  &  C  \\
D  &  E
\end{pmatrix}$ of size $n\times n$
with nonsingular $k\times k$ {\em pivot block} $B=A^{(k)}$, 
 define 
$S=S(A^{(k)},A)=E-DB^{-1}C$,

the {\em Schur complement} of $A^{(k)}$ in $A$,
and
the block factorizations,

\begin{equation}\label{eqgenp}
\begin{aligned}
A=\begin{pmatrix}
I_k  &  O_{k,r}  \\
DB^{-1}  & I_r
\end{pmatrix}
\begin{pmatrix}
B  &  O_{k,r} \\
O_{r,k}  &  S
\end{pmatrix}
\begin{pmatrix}
I_k  &  B^{-1}C  \\
O_{k,r}  & I_r
\end{pmatrix}
\end{aligned}
\end{equation} 
and
\begin{equation}\label{eqgenpin}
\begin{aligned}
A^{-1}=\begin{pmatrix}
I_k  &  -B^{-1}C  \\
O_{k,r}  & I_r
\end{pmatrix}
\begin{pmatrix}
B^{-1}  &  O_{k,r} \\
O_{r,k}  &  S^{-1}
\end{pmatrix}
\begin{pmatrix}
I_k  &  O_{k,r}  \\
-DB^{-1}  & I_r
\end{pmatrix}
\end{aligned}.
\end{equation} 

We verify readily that $S^{-1}$ is the $(n-k)\times (n-k)$ trailing
(that is, southeastern) block of the inverse matrix $A^{-1}$, and so  
the Schur complement $S$ is nonsingular since the matrix $A$ is nonsingular.

Factorization (\ref{eqgenpin}) reduces the inversion of the matrix $A$ 
to the inversion of the leading block $B$ and its 
Schur complement $S$, and we can 
recursively reduce the task to the case of the leading blocks 
and Schur complements of decreasing sizes 
as long as the leading blocks are nonsingular.  
After sufficiently many  recursive steps 
of this
process of
block Gaussian elimination,
we only need to invert matrices  of
small sizes, and then we can stop the process   
 and  apply a selected black box
inversion algorithm. 

In $\lceil\log_2(n)\rceil$ recursive steps
all pivot blocks and 
all other
matrices involved into the 
resulting factorization
turn into scalars, 
all matrix 
multiplications and inversions turn into 
scalar multiplications and divisions,
and we arrive at a
{\em complete recursive factorization} of the matrix $A$.
If
$k=1$ at all recursive steps, then the complete
recursive factorization (\ref{eqgenpin})
defines GENP and can be applied to computing the inverse $A^{-1}$
or the solution ${\bf y}=A^{-1}{\bf b}$
to a linear system $A{\bf y}={\bf b}$. 

Actually, however, any complete recursive factorizations
turns into GENP up to the order in which we 
consider its steps.
This follows because  at most $n-1$  distinct
Schur complements $S=S(A^{(k)},A)$ 
for $k=1,\dots,n-1$
are involved in all recursive block factorization
processes for $n\times n$ matrices $A$,
and so we arrive at the same Schur complement in a fixed position
via GENP and via any other recursive block factorization (\ref{eqgenp}).
Hence we can interpret factorization step  (\ref{eqgenp})
as the block elimination of the
first $k$
columns of the matrix $A$, 

which produces  the  matrix $S=S(A^{(k)},A)$.
If the dimensions 
$d_1,\dots,d_r$ and $\bar d_1,\dots,\bar d_{\bar r}$ of 
the pivot  blocks in 
two block elimination processes
 sum to the same integer $k$, that is, if 
$k=d_1+\dots+d_r=\bar d_1+\dots+\bar d_{\bar r}$,
then 
both processes produce the same Schur complement $S=S(A^{(k)},A)$.
The following results extend this observation. 


\begin{theorem}\label{thsch}
In the recursive block factorization process based on (\ref{eqgenp}),  
 every diagonal block of every block diagonal factor is either 
a leading block of the input matrix $A$ or the Schur complement $S(A^{(h)},A^{(k)})$
 for some integers $h$ and $k$ such that $0<h<k\le n$ and

$S(A^{(h)},A^{(k)})=(S(A^{(h)},A))^{(h)}$.

\end{theorem}


\begin{corollary}\label{corec}
The recursive block factorization process based on equation (\ref{eqgenp})  
can be completed by involving no singular
pivot blocks (and in particular no pivot elements vanish)
if and only if the input matrix $A$ is strongly nonsingular.
\end{corollary}
\begin{proof}
Combine Theorem \ref{thsch} with the equation 
$\det A=(\det B)\det S$,
implied by (\ref{eqgenp}).
\end{proof}


The following theorem bounds the norms of all
 pivot blocks and their inverses and hence bounds
the condition numbers of the blocks, that is,
  precisely the quantities responsible for 
safe numerical performance of block Gaussian elimination 
and GENP. 
		


\begin{theorem}\label{thnorms} (Cf. \cite[Theorem 5.1]{PQZ13}.)
Assume GENP
or block Gaussian elimination
applied to 
 an
$n\times n$
matrix $A$ and
write $N=||A||$ and $N_-=\max_{j=1}^n ||(A^{(j)})^{-1}||$,
and so $N_-N\ge ||A||~||A^{-1}||\ge 1$. 
Then the absolute values of all pivot elements of GENP 
and the norms of all pivot blocks of 
block Gaussian elimination
do not exceed $N_+=N+N_-N^2$,
while the absolute values of the reciprocals of these 
elements and the norms of the inverses of the blocks do not
exceed $N_-$.
\end{theorem}
\begin{proof}
Observe that the inverse $S^{-1}$ of the Schur complement $S$
in (\ref{eqgenp}) is the southeastern block of the inverse $A^{-1}$
and obtain

$||B||\le N$, $||B^{-1}||\le N_-$, and 
$||S^{-1}||\le ||A^{-1}||\le N_-$. Moreover
$||S||\le N+N_-N^2$, due to (\ref{eqgenp}). Now the claimed bound follows 
from Theorem \ref{thsch}.
\end{proof}
\noindent
 


\begin{remark}\label{reir}
By virtue of Theorem \ref{thnorms} 
the norms of the inverses of
all pivot blocks involved into a complete
(and hence also into any incomplete) recursive factorization
of a strongly nonsingular matrix $A$
are at most $N_-$.
We have a reasonable upper bound on $N_-$ if
the matrix $A$ is strongly well-con\-di\-tioned as well.
Then
in view of Theorem \ref{thpert}
the inversion of all pivot blocks is numerically
safe, and we say that {\em GENP is locally safe} for the matrix $A$.


\end{remark}


\begin{remark}\label{renrm}
In the recursive factorizations above 
only the factors of the leading blocks  
and the Schur complements
can contribute to the magnification
of any input perturbation. Namely at most $\lceil\log_2(n)\rceil$
such factors can contribute to the norm of each 
of the output triangular or block triangular
factors $L$ and $U$. This implies the moderately large 
worst case upper bound $(N_+N_-)^{\log_2(n)}$
on their norms, which  
 is overly pessimistic
according to our tests.


\end{remark}


\begin{remark}\label{rerect}
Our study in this and the next two sections can be extended readily to 
the cases of GENP and block Gaussian 
elimination applied to rectangular and possibly rank deficient matrices 
and to under- and over-determined and possibly rank deficient 
linear systems of equations.
Recursive factorization and elimination 
can be completed and are numerically safe when they are applied to any  
strongly nonsingular and  strongly well-conditioned leading block 
of the input matrix, in particular 
to the input matrix 
itself if it is strongly nonsingular and  strongly well-conditioned.
\end{remark}


\section{Singular values of the  matrix products (deterministic estimates)
and GENP and block Gaussian elimination with preprocessing}\label{smrc}


{\em Preprocessing} $A\rightarrow FAH$ 
for a pair of nonsingular matrices $F$ and $H$, 
one of which can be the identity matrix $I$,
reduces
the inversion of a matrix $A$ to
the inversion of a the product $FAH$,
and similarly for the solution of a linear system of equations.  
\begin{fact}\label{faprepr}

Assume three nonsingular matrices $F$, $A$, and $H$
 and a vector ${\bf b}$.
Then 

$A^{-1}=H(AH)^{-1}$, $A^{-1}=(FA)^{-1}F$,
$A^{-1}=H(FAH)^{-1}F$. Moreover,  if $A{\bf x}={\bf b}$, then $AH{\bf y}={\bf b}$, $FA{\bf x}=F{\bf b}$, and 
$FAH{\bf y}=F{\bf b}$, ${\bf x}=H{\bf y}$. 

\end{fact}

 Remark \ref{reir} motivates the choice of the multipliers $F$ 
and $H$ for which the 
matrix $FAH$ is strongly nonsingular and
 strongly well-con\-di\-tioned. 
This is likely to occur already if one of the multipliers
$F$ and $H$
 is the identity matrix 
 and another one is a Gaussian random matrix. 
The studies of 
pre-multiplication by $F$ and post-multiplication by $H$
are similar, and so we only prove the latter claim
in the case of post-multiplication.
We complete our proof in Section \ref{ssgnp}. It
 involves the norms of the inverses of the matrices 
 $(AH)_{k,k}=A_{k,n}H_{n,k}$ for $k=1,\dots,r$,
 which we  estimate in this section  
assuming  
nonrandom multipliers
$H$. 


We begin with  two simple lemmas.


\begin{lemma}\label{lepr2} 
If $S$ and $T$ 
are square orthogonal matrices, then 
$\sigma_{j}(SA)=\sigma_j(AT)=\sigma_j(A)$ for all $j$.
\end{lemma}


\begin{lemma}\label{lepr1} 
Suppose $\Sigma=\diag(\sigma_i)_{i=1}^{n}$, $\sigma_1\ge \sigma_2\ge \cdots \ge \sigma_n$,
 and  $H\in \mathbb R^{n\times r}$.
Then 
$$\sigma_{j}(\Sigma H)\ge\sigma_{j}(H)\sigma_n {\rm for~all}~j.$$
If also $\sigma_n>0$, then 
 $$\rank (\Sigma H)=\rank (H).$$
\end{lemma}


We also need the following basic results (cf. \cite[Corollary 8.6.3]{GL13}).
 
\begin{theorem}\label{thcndsb} 
If $A_0$ is a 
submatrix of a 
matrix $A$, 
then
$\sigma_{j} (A)\ge \sigma_{j} (A_0)$ for all $j$.
\end{theorem} 


 


 


\begin{theorem}\label{thintpr}  
Suppose $r+l\le n\le m$, $l\ge 0$,
 $A\in \mathbb R^{m\times n}$,
$\rank(A_{m,r})=r$ and $\rank(A_{m,r+l})=r+l$.

Then $||A_{m,r}^+||\le ||A_{m,r+l}^+||$.
\end{theorem}


The following theorem will enable us to
estimate the norm $||(AH)^+||$.

\begin{theorem}\label{1}

Suppose $A\in \mathbb R^{n\times n}$,
 $H\in\mathbb R^{n\times r}$,
$\rank (A)=n\ge r$,  
 $A=S_A\Sigma_AT^T_A$ is SVD (cf. (\ref{eqsvd})), 
 and
$\widehat H=T_A^TH$.
 
 
 Then 
\begin{equation}\label{eqah}
\sigma_j(AH)\ge \sigma_l(A) ~\sigma_j(\widehat H_{l,r})~{\rm for~all}~l\le n~{\rm and~all}~j.
\end{equation}
\end{theorem}


\begin{proof}
Note that
 $AH=S_A\Sigma_AT_A^TH$, and so  
$\sigma_j(AH)=\sigma_j(\Sigma_AT_A^TH)=
\sigma_j(\Sigma_A\widehat H)$ for all $j$ 
by virtue of Lemma \ref{lepr2}, because 
$S_A$ is a square orthogonal matrix. 
Moreover it follows from Theorem \ref{thcndsb} that
$\sigma_j(\Sigma_A\widehat H)\ge 
\sigma_j(\Sigma_{l,A}\widehat H_{l,r})$
for all $l\le n$. Combine this bound with the latter equations
and apply Lemma \ref{lepr1}. 
\end{proof}


\begin{corollary}\label{cor1}
Keep the assumptions 

of Theorem \ref{1}. Then 

(i)  $\sigma_{r}(AH)\ge 
\sigma_{\rho}(A)\sigma_{r}(\widehat H_{n,r})=
\sigma_{r}(\widehat H_{n,r})/||A^+||$,

(ii) $||(AH)^+||\le ||A^+||~||\widehat H_{n,r}^+||$
if $\rank(AH)=\rank(\widehat H_{n,r})=r$.
\end{corollary}


\begin{proof}
Substitute $j=r$ and $l=n$
into bound (\ref{eqah}), recall
(\ref{eqnrm+}), and
obtain part (i).
If $\rank(AH)=\rank(\widehat H_{l,r})=r$,
then apply (\ref{eqnrm+})
to obtain that 
$\sigma_{r}(AH)=1/||(AH)^+||$ and 
$\sigma_{r}(\widehat H_{l,r})=1/||\widehat H_{l,r}^+||$. 
Substitute these equations 
into part (i) and obtain part (ii).
\end{proof}


Let us  extend the estimates of Theorem \ref{1}  to the leading
 blocks of a matrix product.


\begin{corollary}\label{cogh}
Keep the assumptions 

of Theorem \ref{1} and also suppose that
 the matrices  $(AH)_{k,k}$ and $\widehat H_{n,k}$ 
have full rank $k$ for
a positive integer
$k\le n$. 

 


Then  
 $$||(AH)_{k,k}^+||\le ||\widehat H_{n,k}^+||~||A_{k,n}^+||\le
||\widehat H_{n,k}^+||~||A^+||.$$ 
\end{corollary} 


\begin{proof}


Note that 
$(AH)_{k,k}=A_{k,n}H_{n,k}$
and that the matrix $A_{k,n}$ has full rank.
Apply Corollary \ref{cor1}  for $A$ and $H$ replaced by 
$A_{k,n}$ and $H_{n,k}$, respectively,
and obtain that $||(AH)_{k,k}^+||\le ||\widehat H_{k,n}^+||~||A_{n,k}^+||$.
Combine (\ref{eqnrm+}) and Theorem \ref{thintpr} and deduce that 
$||A_{n,k}^+||\le||A^+||$. Combine  the two latter 
inequalities 
to complete the proof of part (i). 
Similarly prove part (ii).
\end{proof}


Fact \ref{faprepr}, 
Corollary \ref{corec} and Theorem \ref{thnorms} together imply the
following result.


\begin{corollary}\label{colocsf}
 Suppose that
$A\in \mathbb R^{n\times n}$,
$H\in\mathbb R^{n\times r}$,
 $r\le n=\rank (A)$,  and 
 the matrices $(AH)_{k,k}$ 
 
are   
strongly nonsingular and
   strongly well-con\-di\-tioned
for $k=1,\dots,r$. Then
 GENP  and 
block Gaussian elimination
are locally safe for the
matrix product $AH$
 (see Remark \ref{reir} on the concept ``locally safe").
\end{corollary}


\section{Benefits of using Gaussian multipliers for GENP and block Gaussian elimination


}\label{sapsr1}


 In  Section
\ref{srrm} we recall the norm and condition estimates for  
Gaussian matrices and deduce that these matrices are
strongly nonsingular with probability 1 and
are expected to be strongly well-conditioned.
In Section \ref{ssgnp} we prove
that the pair $(H,A)$
for a nonsingular and well conditioned  matrix $A$ and a
Gaussian matrix $H$ is expected to
satisfy the assumptions of Corollary \ref{colocsf},  
implying that the application of 
GENP  and block Gaussian elimination 
to the product $AH$ is  numerically safe.


\begin{remark}\label{remu0}
The above results do not hold 
if the mean greatly exceeds standard deviation 
of the i.i.d. entries of a
multipliers $H$. Its power for
 achieving numerically safe GENP
and block Gaussian elimination
is usually lost in this case.
Indeed assume 
 a mean $\mu$ and a standard deviation
$\sigma$ such that
$\mu\gg \sigma$ (already   
$\mu>10 \sigma\log (n)$ is suficient).
In this case the matrix $H$ is expected to be 
closely approximated by the rank-1 matrix $\mu {\bf e}{\bf e}^T$
where ${\bf e}^T=(1,1,\dots,1)$.
\end{remark}


\subsection{A Gaussian 

matrix, its rank, norm and condition estimates}\label{srrm}


\begin{definition}\label{defrndm}
A matrix is said to be {\em standard Gaussian random} 
(hereafter we say just
{\em Gaussian}) if it is filled with i.i.d.
Gaussian random
variables  having mean $0$ and variance $1$. 

\end{definition}


\begin{theorem}\label{thdgr1}
A Gaussian matrix $G$ is strongly nonsingular 
with probability 1.
\end{theorem}
\begin{proof}
Assume that the $j\times j$ leading submatrix $G^{(j)}$ of a $k\times l$
Gaussian matrix $G$ is singular for some positive integer $j\le h=\min\{k,l\}$,
that is, $\det(G^{(j)})=0$. 
Since $\det(G^{(j)})$ is a polynomial in the entries of the Gaussian matrix 
 $G^{(j)}$, such matrices form 

an algebraic variety of a lower dimension  in the linear space 
$\mathbb R^{j^2}$. 
($V$ is an algebraic variety of a dimension $d\le N$
in the space $\mathbb  R^{N}$ if it is defined by $N-d$
polynomial equations and cannot be defined by fewer equations.)
Clearly, Lebesgue (uniform) and
Gaussian measures of such a variety equal 0,  being absolutely continuous
with respect to one another. 
Hence these measures of the union of $h$ such matrices

are also 0.
\end{proof}


\begin{theorem}\label{thdgr}
Assume a nonsingular $n\times n$ matrix $A$ and an $n\times k$  
 Gaussian matrix $H_{n,k}$.
Then the product $A_{k,n}H_{n,k}$ is nonsingular with probability 1.
\end{theorem}
\begin{proof}
$\det(A_{k,n}H_{n,k})$ is a polynomials in the entries of the Gaussian matrix $H_{n,k}$.
Such a polynomial vanishes with probability 0 unless it vanishes identically in  $H_{n,k}$,
but  the matrix 
$A_{k,n}A_{k,n}^T$ is positive definite, and so
$\det(A_{k,n}H_{n,k})>0$ for $H_{n,k}=A_{k,n}^T$.
\end{proof}


\begin{definition}\label{defnu}
 $\nu_{j,m,n}$ denotes the random variables
$\sigma_j(G)$ for a Gaussian $m\times n$ matrix $G$
and all $j$, while
$\nu_{m,n}$, $\nu_{F,m,n}$,  $\nu_{m,n}^+$, and $\kappa_{m,n}$ denote the 
 random variables
 $||G||$,  $||G||_F$, 
 $||G^+||$, and $\kappa(G)=||G||~||G^+||$, respectively.
\end{definition}

Note that $\nu_{j,n,m}=\nu_{j,m,n}$, $\nu_{n,m}=\nu_{m,n}$,
$\nu_{n,m}^+=\nu_{m,n}^+$, and $\kappa_{n,m}=\kappa_{m,n}$.

\begin{theorem}\label{thsignorm}
(Cf. \cite[Theorem II.7]{DS01}.)
Suppose 

$h=\max\{m,n\}$, $t\ge 0$,  and
$z\ge 2\sqrt {h}$. 
Then  ${\rm Probability}~\{\nu_{m,n}>z\}\le
\exp(-(z-2\sqrt {h})^2/2\}$ and 
${\rm Probability}~\{\nu_{m,n}>t+\sqrt m+\sqrt n\}\le
\exp(-t^2/2)$.
\end{theorem}


\begin{theorem}\label{thsiguna} 
(Cf. \cite[Proof of Lemma 4.1]{CD05}.)
Suppose 
$m\ge n\ge 2$, and $x>0$ and write $\Gamma(x)=
\int_0^{\infty}\exp(-t)t^{x-1}dt$ 
and $\zeta(t)=

t^{m-1}m^{m/2}2^{(2-m)/2}\exp(-mt^2/2)/\Gamma(m/2)$. 
Then 
 ${\rm Probability}~\{\nu_{m,n}^+\ge m/x^2\}<\frac{x^{m-n+1}}{\Gamma(m-n+2)}$.
\end{theorem}


The following condition estimates
 from \cite[Theorem 4.5]{CD05}

are quite tight for large values $x$, but for $n\ge 2$ even tighter
estimates (although more involved)
 can be found in \cite{ES05}. (See \cite{D88} and  \cite{E88}
on the early study.)


\begin{theorem}\label{thmsiguna}   

If $m\ge n\ge 2$, then 
$${\rm Probability}~\{\kappa_{m,n}m/(m-n+1)>x\}\le \frac{1}{2\pi}(6.414/x)^{m-n+1}$$
for $x\ge m-n+1$, while $\kappa_{m,1}=1$
with probability 1.
\end{theorem}


 


\begin{corollary}\label{cogfrwc} 
 A Gaussian 

matrix  is expected to be  strongly well-con\-di\-tioned.
\end{corollary}


\subsection{Supporting  GENP with Gaussian multipliers


}\label{ssgnp}


The main result of this section  is
Corollary \ref{cogmforalg}, which supports application of 
GENP  and block Gaussian elimination 
to the product $AH$ of a nonsingular matrix $A$ and  a
Gaussian matrix $H$.


We need the following simple basic lemma.


\begin{lemma}\label{lepr3}

Suppose $H$ is a Gaussian matrix, 
$S$ and $T$ are orthogonal matrices, $H\in \mathbb R^{m\times n}$,
$S\in \mathbb R^{k\times m}$, and $T\in \mathbb R^{n\times k}$
for some $k$, $m$, and $n$.
Then $SH$ and $HT$ are Gaussian matrices.
\end{lemma}


\begin{corollary}\label{cor10}
Suppose $A$ is a nonsingular $n\times n$ matrix,
$H$ is an $n\times k$ Gaussian matrix,
for $0<k\le n$,
and $\nu_{g,h}$ and $\nu_{g,h}^+$ are the random values 
of Definition \ref{defnu}. 
Then 

(i) the matrix $(AH)_{k,k})$ is nonsingular 
with probability $1$,

(ii) $||(AH)_{k,k}||\le \nu_{n,k}||A_{k,n}||\le \nu_{n,k}||A||$, and

 (iii) $||(AH)_{k,k}^{+}||\leq \nu_{n,k}^+||A^{+}||$.

\end{corollary}


\begin{proof}
Part (i) restates Theorem \ref{thdgr}.
Part (ii) follows because $(AH)_{k,k}=A_{k,n}H_{n,k}$,
 $H_{n,k}$ is a Gaussian  matrix,
 and $||A_{k,n}||\le ||A||$. 
Part (iii) follows from Corollary \ref{cogh} 
because 
 $\widehat H_{n,k}$ is a Gaussian  matrix
by virtue of Lemma \ref{lepr3}. 
\end{proof}


\begin{corollary}\label{cogmforalg}
Suppose that the $n\times n$ matrix $A$
is nonsingular and  well-con\-di\-tioned.
Then the  choice of 
Gaussian multiplier 
$H$  is expected 
to satisfy the assumptions of Corollary \ref{colocsf}.
\end{corollary}
\begin{proof}
Recall that $(AH)_{k,k}=A_{k,n}H_{n,k}$
and hence $||(AH)_{k,k}||=||A_{k,n}||\nu_{n,k}$.
Then
combine Theorems 
\ref{thdgr},
\ref{thsignorm}, and
\ref{thsiguna} 
and Corollary
 \ref{cor10}.
\end{proof}


\section{Random structured multipliers for GENP and block Gaussian elimination}\label{srnd}


This subsection involves complex  matrices. A complex matrix $M$ is  unitary if
 $M^HM=I$ or $MM^H=I$ where $M^H$ denotes its Hermitian transpose,
so that $M^H=M^T$ for a real matrix $M$. 

Hereafter 
 $\omega=\omega_n=\exp(\frac{2\pi}{n} \sqrt {-1})$ 
denotes an $n$th primitive root of unity,
$\Omega=(\omega^{ij})_{i,j=0}^{n-1}$ 
is the matrix of the discrete Fourier transform 
 at $n$ points
(we use the acronym {\em DFT}),
and $\Omega^{-1}=\frac{1}{n}\Omega^{H}$. 

An $n\times n$
circulant matrix $C=(c_{i-j\mod n})_{i,j=0}^{n-1}$ 
is defined by its first column ${\bf c}=(c_i)_{i=0}^{n-1}$.

\begin{example}\label{exc1} {\em  Generation of random real circulant
  matrices.}
Generate the vector ${\bf c}$ of $n$ i.i.d. random real variables 
in the range $[-1,1]$ under the uniform probability distribution
on this range.
Define an $n\times n$ circulant matrix $C$ with the first column ${\bf c}$.
\end{example}

The following theorem links 
the matrices $\Omega$ and $\Omega^{-1}$ to the  class of circulant matrices.

\begin{theorem}\label{thcpw} (Cf. \cite{CPW74}.)
 Let
$C$ denote a circulant $n\times n$ matrix defined by its
first column ${\bf c}$ and write 
${\bf u}=(u_i)_{i=1}^n=\Omega {\bf c}$.
Then
$C=\Omega^{-1}\diag(u_j)_{j=1}^n\Omega$. Furthermore 
$C^{-1}=\Omega^{-1}\diag(1/u_j)_{j=1}^n\Omega$ if 
the matrix $C$ is nonsingular.
\end{theorem}
By using FFT, one can multiply the matrices $\Omega$
and $\Omega^H=\Omega^{-1}$ by a vector by using 
$O(n\log (n))$ flops for any $n$ (cf., e.g., \cite[page 29]{p01}), and 
Theorem \ref{thcpw} extends this complexity bound to 
multiplication  of an
$n\times n$ 
circulant matrix and its inverses by a vector.

We need $2n^3-n^2$ flops in order to compute the product $AH$  of
the pair of $n\times n$ matrices $A$  
and  $H$.
If, however, $H$ is a   
 circulant matrix, then 
we can compute $AH$ by using order of $n^2\log (n)$ flops.
For a Toeplitz-like matrix $A$
defined by its displacement generator of bounded length $l$,
we use  $O(ln\log(n))$ flops in order to compute a
displacement generator of length $l$ for the matrix $AH$.
(See \cite{p01} for the definition of displacement generators.)
In the case of Toeplitz matrices we have $l\le 2$ and use
 $O(n\log(n))$ flops.
This motivates using Gaussian circulant multipliers $H$,
that is, circulant matrices $H$ whose first column vector is 
Gaussian. It has been proved in \cite{PSZa}
that such matrices are expected to be well-conditioned,
which is required for any multiplicative preconditioner.


We can define a unitary circulant matrix
by its first column vector
${\bf c}=\Omega (\exp(r_i\sqrt {-1}))_{i=0}^{n-1}$
for any set of real values $r_0,\dots,r_{n-1}$.

\begin{example}\label{exunc} {\em Generation of  random unitary circulant 
 matrices.}

(i) Generate a vector ${\bf u}=(u_j)_{j=1}^n$ where
 $u_j=\exp(2\pi \phi_j \sqrt {-1})$ (and so
$|u_j|=1$ for all $i$) and  where $\phi_1,\dots,\phi_n$
 are $n$ independent random real variables, e.g.,
Gaussian  variables or the variables uniformly distributed in the range $[0,1)$. 

(ii) Compute the vector  
${\bf c}=\Omega^{-1}{\bf u}$,
where  $\Omega$
denotes the $n\times n$ DFT  matrix.
Output 
the unitary circulant matrix $C$ 
defined by its first column ${\bf c}$.
\end{example}

Our proof that Gaussian multipliers enforce strong nonsingularity 
of a nonsingular matrix with probability 1 (see Theorem \ref{thdgr}) 
has been non-trivially extended in \cite{PZa} to the case of 
Gaussian circulant multipliers.
Furthermore 
strong nonsingularity  holds
with probability close to 1 if we 
fill the first column of a multiplier $F$ or $H$ 
with i.i.d. random variables defined
under the uniform probability distribution
over a sufficiently large finite set
(see Appendix \ref{srsnrm} and  \cite{PSZa}).

In our tests with random input matrices, 
Gaussian circulant and general Gaussian multipliers have shown 
the same power  
of supporting numerically safe GENP 
(see Section \ref{sexgeneral}), but we cannot extend
our basic Lemma \ref{lepr3} and
our Corollary \ref{cogmforalg}
to the  case of circulant matrices.
Moreover our Theorem \ref{thcgenp} and
Remark \ref{remext}  below
show that, for a  specific narrow class of input matrices $A$, 
GENP with these multipliers is expected to fail numerically. 

\begin{theorem}\label{thcgenp}
Assume a large integer $n$ and the $n\times n$ DFT matrix $\Omega$,
which is unitary up to scaling by $1/\sqrt n$. 

(i) Then application of GENP to this matrix fails numerically and
 
(ii) a Gaussian circulant $n\times n$ multiplier $C=\Omega^{-1}D\Omega$
with Gaussian diagonal matrix $D=\diag(g_j)_{j=1}^n$
(having i.i.d. Gaussian diagonal entries $g_1,\dots,g_n$)
is not expected to fix this problem. 
\end{theorem}
\begin{proof}
 (i) Subtract  the 
first row of the block 
$\Omega_{2,2}$ of the matrix $\Omega$ 
and  the resulting vector
from its second row. Obtain the vector
 $(0,\omega-1)$ with the norm $|\omega-1|=2\sin (\pi/n)$.
Assume that $n$ is large and then observe that 
$2\sin (\pi/n)\approx 2\pi/n$ and that
 the  variable $2\pi/n$ is expected to be small, implying that
 $\nrank(\Omega_{2,2})=1$ because $||\Omega_{2,2}||\ge \sqrt 2$
for large $n$. 

(ii) Note that $\Omega C=D \Omega$.
The Gaussian variable $g_1$
vanishes with probability 0, 
and so we can assume that $g_1\neq 0$.
 Multiply  the 
first row of the block 
$(D\Omega)_{2,2}$ of the matrix $D\Omega$ by $g_2/g_1$
and subtract the resulting vector
from the second row. Obtain the vector
 $(0,(\omega-1)g_2)$ with the norm
 $|(\omega-1)g_2|=2|g_2\sin (\pi/n)|$
Assume that $n$ is large and then observe that 
$|(\omega-1)g_2|\approx 2|g_2|\pi/n$
and that
 the  variable $2|g_2|\pi/n$ is expected to be small. 
Hence
 $\nrank((\Omega C)_{2,2})=\nrank((D\Omega)_{2,2})$ is expected to equal 1
 because 
$||\Omega C)_{2,2}||\le ||\Omega_{2,2}||\max\{g_1,g_2\}$,
$||\Omega_{2,2}||\ge \sqrt 2$
 and the random variable  $\max\{|g_1|,|g_2|\}$
is not expected to be close to 0.
\end{proof}

\begin{remark}\label{remext}
The same argument shows that 
Gaussian circulant multipliers $C$ are not expected to support GENP
for a bit larger class of matrices, e.g., for $A=M\Omega$ where $M=\diag(D_i)_{i=1}^k$,
 $D_1=\diag(d_1,d_2)$,  and $d_1$  and $d_2$ are two positive constants
and the input size $n\times n$ is large 
as well as where the matrix  $M$ is strongly diagonally dominant. 
The reader is challenged to 
find out whether   
GENP with a Gaussian circulant preprocessor is 
expected to fail numerically for other  classes of input matrices,
in particular for any
  subclass of the classes of Toeplitz or
Toeplitz-like matrices (cf. \cite{p01} and \cite{P15} on these classes).
Another challenge is to choose a distinct random structured preprocessor
for which the above problem is avoided. E.g., consider
the product $\prod_{i=1}^h C_i$ where $h$ is 
 a small integer exceeding 1,
$C_{2j}$
are circulant matrices and $C_{2j-1}$ are skew-circulant 
(see the definition in \cite{p01}).
 Toward the same goal
 we can apply simultaneously random structured
 pre- and post-multipliers $F$ and $H$,
defined by some i.i.d. random parameters, or the pairs of 
PRMB multipliers of \cite{BBD12}. In the case of Toeplitz or Toeplitz-like
input matrices $A$, the multiplications $FA$ and $AH$ are much 
less costly if the multipliers $F$ and $H$ are
circulant matrices, skew-circulant matrices, or the products of 
such matrices.
\end{remark}


\section{Low-rank approximation


}\label{sapsr}


Suppose we seek a rank-$r$ approximation of a matrix $ A$ 
that has a small numerical rank $r$.
One can solve this problem by computing 
 SVD of the 
 matrix $A$
 or, at a lower cost, by computing
 its rank-revealing 
 factorization \cite{GE96}, \cite{HP92}, \cite{P00a},
but using
random matrix multipliers instead
has some benefits \cite{HMT11}.
In this section we study the latter randomized approach. In
 its first subsection we recall some relevant definitions and 
auxiliary results. 


\subsection{Truncation of SVD. Leading and trailing singular spaces}\label{sosvdi1}


Truncate the square orthogonal matrices $S_A$ and $T_A$
and the square diagonal matrix $\Sigma_{A}$ of the SVD of (\ref{eqsvd}),
write $S_{\rho,A}=(S_A)_{m,\rho}$,
 $T_{\rho,A}=(T_A)_{n,\rho}$, and  $\Sigma_{\rho,A}=(\Sigma_{A})_{\rho,\rho}=\diag(\sigma_j)_{j=1}^{\rho}$,
and obtain {\em thin SVD}
\begin{equation}\label{eqthn}
A=S_{\rho,A}\Sigma_{\rho,A}T_{\rho,A}^T,~\rho=\rank (A).
\end{equation}
Now for every integer $r$ in the range $1\le r\le\rho=\rank(A)$,  
write $\Sigma_{\rho,A}=\diag(\Sigma_{r,A},\bar\Sigma_{A,r})$ and
partition
the matrices $S_{\rho,A}$ and $T_{\rho,A}$ into block columns,
 $S_{\rho,A}=(S_{r,A}~|~\bar S_{A,r})$, and 
 $T_{\rho,A}=(T_{r,A}~|~\bar T_{A,r})$
where $\Sigma_{r,A}=(\Sigma_{A})_{r,r}=\diag(\sigma_j)_{j=1}^r$,
$S_{r,A}=(S_A)_{m,r}$,
and $T_{r,A}=(T_A)_{n,r}$.
Then  
 partition the thin SVD as follows,
\begin{equation}\label{eqldtr}
A_r=S_{r,A}\Sigma_{r,A}T_{r,A}^T,~\bar A_r=\bar S_{A,r}\bar \Sigma_{A,r}\bar T_{A,r}^T,~
A=A_r+\bar A_r~{\rm for}~ 1\le r\le\rho=\rank(A),
\end{equation} 
and call the above decomposition 

the $r$-{\em truncation of  thin SVD} (\ref{eqthn}).
Note that $\bar A_{\rho}$ is an empty matrix
and recall that
\begin{equation}\label{eqsvdpert} 
||A-A_r||= \sigma_{r+1}(A).
\end{equation}
Let
$\mathbb S_{r,A}$ and  $\mathbb T_{r,A}$ denote
the ranges (that is, the column spans)
 of the matrices $S_{r,A}$ and $T_{A,r}$, respectively.
If $\sigma_r>\sigma_{r+1}$, 
then 
$\mathbb S_{r,A}$ and $\mathbb T_{r,A}$ are
the left and right {\em leading singular spaces}, respectively,
    associated with the $r$ largest singular values of the matrix $A$.


The left singular spaces of a matrix $A$ are 
the right   singular spaces of its transpose $A^T$ and vice versa.
All matrix bases for the singular spaces $\mathbb S_{r,A}$ and $\mathbb T_{r,A}$
are given by the matrices $S_{r,A}X$ and
$T_{r,A}Y$, respectively,
for  nonsingular $r\times r$ matrices $X$ and $Y$.
The bases are 
orthogonal if the matrices $X$ and $Y$ 
are orthogonal.


\subsection{The basic algorithm}\label{sprb0}


Assume an $m\times n$ matrix $A$ 
having a small numerical rank $r$ and 
 a Gaussian $n\times r$
matrix $H$. Then according to 
\cite[Theorem 4.1]{HMT11}, 
the column span of the matrices $AH$ and  
$Q(AH)$ is likely to approximate the leading singular space $\mathbb S_{r,A}$
of the matrix $A$,  
and if it does, then it follows that the rank-$r$
matrix $QQ^TA$ approximates the matrix $A$.

In this subsection we recall the algorithm
supporting this theorem,  where temporarily we assume
nonrandom  multipliers $H$.
 In the next subsections we keep it nonrandom and
estimate the output approximation errors of the algorithm assuming
no oversampling,
suggested in \cite{HMT11}.  Then we
extend our study to the case where $H$ is a Gaussian, 
and  in Section \ref{sapgm} cover 
the results in the case of random structured 
multipliers.


\begin{algorithm}\label{algbasmp} {\bf Low-rank approximation of a matrix.} 
(Cf. Remarks \ref{reaccr} and \ref{relss}.)
 


\begin{description}


\item[{\sc Input:}] 
A
  matrix 
$ A\in \mathbb R^{m\times n}$,
its 
numerical rank $r$, and
two integers $p\ge 2$ and $l=r+p\le \min\{m,n\}$.


\item[{\sc Output:}]
an orthogonal matrix $Q\in \mathbb R^{m\times l}$ such that

the matrix
$QQ^TA\in \mathbb R^{m\times n}$ has rank 
at most $l$ and approximates the matrix $A$.


\item[{\sc Initialization:}] $~$
Generate  an
$n\times l$
 matrix $H$.


\item[{\sc Computations:}] $~$


\begin{enumerate}


\item 
Compute an $n\times l$ orthogonal 
matrix $Q=Q(AH)$, sharing its range with 
the matrix $AH$.

\item 
Compute


and output the matrix $R_{AH}A=QQ^TA$ and stop. 


\end{enumerate}


\end{description}


\end{algorithm}

This basic algorithm from \cite{HMT11} uses $O(lmn)$ flops overall.
 


\subsection{Analysis of the basic algorithm assuming 
no randomization and no oversampling}\label{saaba}


In Corollaries \ref{coerrb} and  \ref{coover} of this subsection 
we estimate the error norms for the approximations
computed by Algorithm \ref{algbasmp} whose
 oversampling parameter $p$ is set to 0,
namely for the  
approximation of an orthogonal  basis for the leading singular space
$\mathbb S_{r,A}$ (by column set of the matrix $Q$ of the algorithm)
and for a
rank-$r$ approximation of the matrix $A$.
We first recall the following results.

\begin{theorem}\label{thsng} (Cf.  (\ref{eqopr}).)
Suppose  $A$ is
an $m\times n$ matrix,
$S_A\Sigma_AT_A^T$ is its SVD, 

 $r$ is an integer,
$0<r\le l\le\min\{m,n\}$, and 
$Q=Q_{r,A}$ is an orthogonal 

matrix 
basis for the space 
$\mathbb S_{r,A}$.
Then 

$||A-QQ^TA||=\sigma_{r+1}(A)$.

\end{theorem}


\begin{theorem}\label{thover} 
Assume two


 
matrices $A\in \mathbb R^{m\times n}$
and $H\in \mathbb R^{n\times r}$ and define the two matrices 
$A_r$ and $\bar A_r$ of (\ref{eqldtr}). Then 
 $AH=A_rH+\bar A_rH$ where $A_rH=S_{r,A}U$, $U=\Sigma_{r,A}T_{r,A}^TH$.
Furthermore
 the columns of the matrix $A_rH$ 
span the space $\mathbb S_{r,A}$
if $\rank (A_rH)=r$.
\end{theorem}


These results together imply that 
the columns of the matrix $Q(AH)$ form an approximate 
orthogonal basis
of the linear space $\mathbb S_A$, and next we
estimate the error norms of this approximations.

\begin{theorem}\label{thuvwz} 
Keep the assumptions of Theorem \ref{thover}.  
Then

(i) $||\bar A_rH||_F\le \sigma_{r+1}(A)~||H||_F$.  

(ii) Furthermore if the matrix $T_{r,A}^TH$ is nonsingular, then
 $||(A_rH)^+||\le ||(T_{r,A}^TH)^{-1}||/\sigma_{r}(A)$. 
\end{theorem}


\begin{proof}
Recall that
\begin{equation}\label{eqfrobun}
||U||=||U||_F=1,~||UAV||\le||A||,
~{\rm and}~||UAV||_F\le||A||_F~{\rm for~orthogonal~matrices}~U~{\rm and}~V. 
\end{equation} 
Then note that  
$||\bar A_rH||_F= ||\bar S_{A,r}\bar \Sigma_{A,r}\bar T_{A,r}^TH||_F\le
||\bar \Sigma_{A,r}\bar T_{A,r}^TH||_F$ by virtue of  bound (\ref{eqfrobun}).

Combine this bound with 
 Lemma \ref{lepr1}  and obtain that 
$||\bar A_rH||_F\le \sigma_{r+1}(A)~||\bar T_{A,r}^TH||_F$,
which is not greater than $\sigma_{r+1}(A)~||H||_F$
by virtue of  bound (\ref{eqfrobun}).
This proves part (i). 

Part (ii) follows because 
$(A_rH)^+=(S_{r,A}\Sigma_{r,A} T_{r,A}^TH)^{-1}=
(T_{r,A}^TH)^{-1}\Sigma_{r,A}^{-1}S_{r,A}^T$ if the matrix $T_{r,A}^TH$ is nonsingular
and because $||S_{r,A}||=1$ and $||\Sigma_{r,A}^{-1}||=1/\sigma_{r}(A)$. 
\end{proof}


Combine  Theorems \ref{thpert1}, \ref{thover}, and \ref{thuvwz} to obtain the following estimates.

\begin{corollary}\label{coerrb} 
Keep the assumptions of Theorem \ref{thover}, let the matrix $T_{r,A}^TH$ be nonsingular
and write $$||E||_F=\sigma_{r+1}(A)~||H||_F,$$ 
$$\Delta_+=\sqrt 2~||E||_F~||(T_{r,A}^TH)^{-1}||~/\sigma_{r}(A)= 
\Delta_+=\sqrt 2~||H||_F~||(T_{r,A}^TH)^{-1}||~\sigma_{r+1}(A)/\sigma_{r}(A).$$
Then 
$$\Delta=||Q(A_rH)^T-Q(AH)^T|| \le \Delta_++O(||E||_F^2).$$ 
\end{corollary} 

 Next combine 
 Corollary \ref{copert} with Theorem \ref{thsng} and
employ the orthogonal projection $P_{AH}=Q(AH)Q(AH)^T$ (cf. (\ref{eqopr}))
to extend the latter estimate to bound the error norm of low-rank approximation of a matrix $A$ 
by means of Algorithm \ref{algbasmp}. 


\begin{corollary}\label{coover} 
Keep the assumptions of Corollary \ref{coerrb} and write

$\Delta_+'= \sigma_{r+1}(A)+2\Delta_+ ||A||$.

 Then 

$$\Delta'=||A-P_{AH}A||\le \Delta_+'+O(||E||_F^2||A||).$$

\end{corollary}


\begin{proof}
Note that 

$||A-P_{AH}A||\le ||A-P_{M}A||+||(P_{M}-P_{AH})A||$
for any $m\times r$ matrix $M$.
Write $M=A_rH$, apply Theorem \ref{thsng}

and  obtain
$||A-P_{M}A||=\sigma_{r+1}(A)$.
Corollaries \ref{copert}
and \ref{coerrb} together imply that
$||(P_{M}-P_{AH})A||\le ||A||~||P_{A_{r}H}-P_{AH}||\le 2 \Delta ||A||$.
 Combine the above relationships.

\end{proof}


\begin{remark}\label{reaccr} 
Write $B_i=(A^TA)^iA$ and 
recall that $\sigma_j(B_i)=(\sigma_j(A))^{2i+1}$ for all positive integers
$i$ and $j$. Therefore
one can 
apply the power transforms
$A\rightarrow B_i$ 
for $i=1,2,\dots$ to
increase the
 ratio $\sigma_{r}(A)/\sigma_{r+1}(A)$, which
shows the gap between the two singular values.
Consequently
 the bound $\Delta_+$ on the error norm  of
the approximation of an orthogonal basis 
 of the leading singular space
$\mathbb S_{r,A}$ by $Q(B_iH)$ 
is expected to decrease as $i$
increases  (cf. \cite[equation (4.5)]{HMT11}).
We use the matrix $AH=B_0H$ in
Algorithm \ref{algbasmp}, but suppose we replace it with 
the matrices $B_iH$ for small positive integer $i$, or even for $i=1$,
which would amount just to symmetrization. Then 
we would obtain low-rank approximation with the optimum error
$\sigma_{r+1}(A)$ up to the terms of higher order in 
$\sigma_{r+1}(A)/\sigma_{r}(A)$
as long as the value $||H||_F||(T_{r,A}^TH)^{-1}||$
is reasonably bounded from above.  
The power transform
$A=B_0\rightarrow B_i$ requires to increase by a factor of $2i+1$
the number of ma\-trix-by-vec\-tor
multiplications involved, but for small positive integers $i$,
the additional  computational cost
is still dominated 
by  the costs of computing the SVD and rank-revealing factorizations.

\end{remark}


\begin{remark}\label{relss}
Let us summarize our analysis. Suppose
that
the ratio $\sigma_{r}(A)/\sigma_{r+1}(A)$
is large and that the matrix product
 $P=T_{r,A}^TH$ has full rank $r$
and is 

well-con\-di\-tioned.
Now set to 0
 the oversampling integer parameter  $p$  
of Algorithm \ref{algbasmp}.
Then,
by virtue of 
Theorem \ref{thuvwz} and
Corollaries \ref{coerrb} and \ref{coover},
 the algorithm 
 outputs a close approximation $Q(AH)$
to an orthogonal bases for the leading singular space
$\mathbb S_{r,A}$
of the input matrix $A$ 
and  a rank-$r$ approximation
to this matrix. Up to the terms of higher order, the error norm 
of the latter approximation
is within a factor of $1+||H||_F~||(T_{r,A}^TH)^{-1}||/\sigma_{r}(A)$
from the optimal bound $\sigma_{r+1}(A)$. 
By applying the above power transform
of the input matrix $A$  
at a low computational  cost, 
 we can  decrease the error norm even below the value
 $\sigma_{r+1}(A)$.
\end{remark}


\subsection{Supporting low-rank approximation with Gaussian multipliers}\label{sapg}


In this subsection we extend the results of the previous one to support
the choice of Gaussian multiplier $H$ in Algorithm \ref{algbasmp}, whose 
``actual outcome 

is very close to the typical outcome because of the measure 
concentration effect" \cite[page 226]{HMT11}.

\begin{theorem}\label{thrrk}
Suppose $A\in \mathbb R^{m\times n}$,
$A=S_A\Sigma_AT_A^T$ is its SVD of (\ref{eqsvd}),
$H=\mathbb R^{n\times r}$
is a Gaussian 
matrix, and $\rank(A)=\rho\ge r$.

(i) Then the matrix $T_{r,A}^TH$ is Gaussian.

(ii) Assume the values $\Delta_+$ and $\Delta_+'$ 
of Corollaries \ref{coerrb} and \ref{coover}
and the values $\nu_{F,n,r}$ and $\nu_{r,r}^+$
of Definition \ref{defnu}.
Then $\Delta_+=\sqrt 2~\nu_{F,n,r} \nu_{r,r}^+\sigma_{r+1}(A)/\sigma_{r}(A)~{\rm and}
 ~\Delta_+'= \sigma_{r+1}(A)+2\Delta_+ ||A||$.

\end{theorem}


\begin{proof}
 $T_{A}^TH$ is a Gaussian matrix by virtue of Lemma \ref{lepr3}.
Therefore so is its square submatrix $T_{r,A}^TH$ as well.  This proves part (i),
which implies part (ii). 
\end{proof}


\begin{corollary}\label{coalg1}
A
Gaussian multiplier 

$H$ is expected 
to 
support safe numerical application of 
  Algorithm \ref{algbasmp} even
 where the oversampling integer parameter $p$ is set to 0.
\end{corollary}
\begin{proof}
Combine Theorems \ref{thdgr1} 
 and \ref{thrrk} with
Corollary \ref{cogfrwc}.
\end{proof}


\subsection{Supporting low-rank approximation with random structured multipliers}\label{sapgm}


 Multiplication of an $n\times n$ matrix $A$ by a
 Gaussian matrix $H$ at Stage 1
of Algorithm \ref{algbasmp} requires $(2r-1)n^2$ flops, 
but  
we can save a factor of $r/\log (r)$ flops by applying
 structured random multipliers $H$.
In particular we can use 
 subsampled random Fourier transforms (SRFTs) of
\cite[equation (4.6)]{HMT11},
subsampled random Hadamard transforms (SRHTs) of
 \cite{T11}, 
the chains of Givens rotations (CGRs) of
 \cite[Remark 4.5.1]{HMT11}, and the leading Toeplitz submatrices
 $C_{n,r}$ and $C_{r,n}$ of 
random circulant $n\times n$ matrices $C$.  
We need just $n$ random parameters to define 
a Gaussian circulant  $n\times n$ matrix $C$
and its leading Toeplitz blocks $C_{n,r}$ and $C_{r,n}$, 
and similarly for the other listed classes of structured matrices.


\begin{example}\label{exsrft}
For two fixed integers $l$ and $n$, $1<l<n$, 
 SRFT  $n\times l$ matrices
are the matrices of the form $S=\sqrt{n/l}~D\Omega R$. Here 
$D$ is a random $n\times n$ diagonal
matrix whose diagonal entries are i.i.d. variables uniformly distributed on the 
unit circle $C(0,1)=\{x:~|x|=1\}$,  $\Omega$ is the DFT matrix,
and $R$ is a random $n\times l$
permutation matrix defined by random choice of $l$ columns under the uniform probability distribution
on the set of the $n$ columns
of the identity matrix $I_n$ (cf. \cite[equation (4.6) and Section 11]{HMT11}).
\end{example}

 Theorem \ref{thcpw} implies the following fact.

\begin{corollary}\label{cocpw}
Assume an $n\times l$ SRFT matrix $S$. Then 
 $\sqrt{l/n}~\Omega^{-1}S$ is an $n\times l$ submatrix of
a unitary circulant $n\times n$ matrix.
\end{corollary}
 
According to the extensive tests by many researchers,
various random structured $n\times l$ multipliers (such as SRFT, SRHT, CGR and CHR matrices)
support low-rank approximation already where 
the oversampling parameter $p=l-r$ is a reasonable constant
(see \cite{HMT11} and \cite{M11}).
In  particular 
 SRFT with oversampling by 20 is adequate in almost all
applications of low-rank approximations \cite[page 279]{HMT11}.
Likewise, in our extensive tests 
covered in Section \ref{stails}, 
Toeplitz multipliers defined as the $n\times r$ leading blocks of
$n\times n$ random circulant 
matrices
consistently  
 supported low-rank approximation
 without oversampling
as efficiently as Gaussian multipliers. 

As in the case of  our
randomized support for GENP and block Gaussian elimination,
formal analysis of
the impact of random structured multipliers 
is complicated because we cannot  use  Lemma \ref{lepr3}.
Nevertheless, by allowing substantial  oversampling,
one can still prove that SRFT multipliers 
are expected to support low-rank approximation 
of a matrix having  a small numerical rank.
     
\begin{theorem}\label{thsrtft}  {\rm Error bounds for
low-rank approximation with SRFT} (cf. \cite[Theorem 11.2]{HMT11}). 
Fix four integers $l$, $m$, $n$, and $r$ 
such that $4[\sqrt r+\sqrt {8\log(rn)n}]^2\log (r)\le l\le n$. 
Assume  an $m\times n$ matrix $A$ with singular
values $\sigma_1\ge\sigma_2\ge \sigma_3\ge \dots$, 
an $n\times l$ SRFT 
matrix $S$ of Example \ref{exsrft},
and $Y = AS$. 
Then with a probability  $1-O(1/r)$ it holds that
$$||(I-P_Y )A|| \le \sqrt{1+7n/l}~\sigma_{r+1}~{\rm and}
~||(I-P_Y )A||_F \le \sqrt{1+7n/l}~(\sum_{j>r} ~\sigma^2_j)^{1/2}.$$
\end{theorem}
 
\begin{remark}\label{remsrtft}
Clearly the theorem still holds
if we  replace the matrix $S$ by the matrix $US$ for 
a unitary matrix $U=(1/\sqrt n)\Omega^{-1}$.
In this case $US=CR$ for the matrix $R$ of Example 
\ref{exsrft} and the circulant matrix $C=\Omega^{-1}D\Omega$ 
 (cf. Theorem \ref{thcpw}).
By virtue of Theorem \ref{thsrtft} 
we can expect that Algorithm \ref{algbasmp} would produce a rank-$r$ approximation 
if we choose a multiplier $H$ being 
an SRFT $n\times l$ matrix  
or
the  $n\times l$ submatrix $CP$ of $n\times n$
 random unitary circulant matrix $C$ made up of its $l$ randomly selected
columns where the selection is defined by the matrix $P$
of Example \ref{exsrft} and where $l$ is an 

integer of order $r\log(r)$. 
Recall that multiplication of an $n\times n$ Toeplitz matrix 
by an $n\times l$ matrix $US=CP$ involves $O(nl \log(n))$ flops \cite{p01},
versus $O(n^2l)$ in the straightforward algorithm.
\end{remark}


\section{Numerical Experiments}\label{sexp}


We performed numerical experiments with  random general,  circulant and Toeplitz  
matrices  by using MATLAB in the Graduate Center of the City University of New York 
on a Dell computer with a Intel Core 2 2.50 GHz processor and 4G memory running 
Windows 7. In particular we generated
Gaussian matrices by using the standard normal distribution function randn of MATLAB,
and we use the MATLAB function rand for generating numbers in the range $[0,1]$
under the uniform probability distribution function for Example \ref{exc1}. 

We display our estimates obtained in terms of the spectral matrix norm but 
our tests showed similar results where we used the Frobenius norm  instead.


\subsection{GENP with Gaussian and random circulant multipliers}\label{sexgeneral}


 
We applied both GENP and the
preprocessed GENP to $n\times n$ DFT matrices $A=\Omega$
and to the matrices $A$ generate as follows.
We fixed $n=2^s$ and $k=n/2$ for $s=6,7,8,9,10$, and first,
by following \cite[Section 28.3]{H02},
generated a $k\times k$ matrix $A_k=U\Sigma V^T$ where we chose  $\Sigma=\diag (\sigma_i)_{i=1}^k$ with
$\sigma_i=1$ for $i=1,\dots,k-4$ and $\sigma_i=0$ for $i=k-3,\dots,k$ and 
where $U$ and $V$ were $k\times k$  random orthonormal  
matrices, computed as 
the $k\times k$ factors $Q(X)$ in the QR factorization of $k\times k$ random matrices $X$. 
Then we
 generated Gaussian  Toeplitz  matrices 
$B$, $C$ and $D$  such that 
 $||B||\approx ||C||\approx ||D||\approx ||A_k||\approx 1$
and defined 
the $n\times n$ matrix  
 $A=\begin{pmatrix}
A_k  &  B  \\
C    &  D
\end{pmatrix}.$
For every dimension $n$, $n=64, 128, 256, 512, 1024$  we run 1000 numerical tests
where we solved the linear system $A{\bf x}={\bf b}$ with  
Gaussian vector ${\bf b}$ and 
output the maximum, minimum and average relative residual norms 
$||A{\bf y}-{\bf b}||/||{\bf b}||$ as well as the standard deviation.
Figure \ref{fig:fig1} and Table \ref{tab61} show the norms of $A^{-1}$. 
They ranged from $2.2\times 10^1$ to  
$3.8\times 10^6$ in our tests. 

At first we describe the results of our tests for the latter class of matrices $A$.
As we expected GEPP has always output accurate solutions 
to the linear systems $A{\bf y}={\bf b}$
 in our tests (see Table \ref{tab62}). 
GENP, however, was expected to fail for these systems, 
because the $(n/2)\times (n/2)$   leading principal block $A_k$ of the matrix $A$  
was singular, having nullity 
$k-\rank (A_k)=4$.
Indeed  this caused poor performance of GENP in our tests,
which have consistently output corrupted solutions,   with  
relative 
residual norms
ranging from $10^{-3}$ to $10^2$. 

In view of Corollary \ref{cogmforalg} 
we expected to fix this deficiency by means of multiplication
by Gaussian matrices, and indeed  in all our tests
we observed 
residual norms below
 $1.3\times 10^{-6}$, and they decreased below $3.6\times 10^{-12}$ 
in a single step of iterative refinement
(see Table \ref{tab63}). Furthermore the tests 
showed the same power of preconditioning 
where we used the circulant multipliers of Examples \ref{exc1} and 
\ref{exunc} (see Tables  \ref{tab64} and \ref{tab65}). 
As can be expected, the output accuracy of GENP with 
preprocessing has deteriorated a little versus GEPP
 in our tests.
The output residual norms, however, were small enough to support
application of the inexpensive iterative refinement. Already its
single step decreased the average relative residual norm 
below $10^{-11}$ for $n=1024$ in all our tests with Gaussian multipliers
and to about $10^{-13}$ for $n=1024$ in all our tests with circulant multipliers
of Examples \ref{exc1} and \ref{exunc}.


 See further details in Figures \ref{fig:fig2} and \ref{fig:fig3} and  Tables \ref{tab63}--\ref{tab65}. This indicates that GENP with preprocessing 
followed by even a single step of iterative refinement 
is backward stable, similarly to the celebrated result of
\cite{S80}.

We also applied similar tests to the
$n\times n$ DFT matrix $A=\Omega$. The results were in very good accordance with our study
in Section \ref{srnd}. Of course in this case the solution of a linear system $A{\bf x}={\bf b}$
can be computed immediately as ${\bf x}=\frac{1}{n}\Omega^H{\bf b}$,
but we were  not seeking the solution, but were trying to
 compare the performance of GENP with and without preprocessing.
In these tests the norm $||A^{-1}||$ was fixed at $1/\sqrt n$.
 GEPP produced the solution within the relative residual norm between 
$10^{-15}$ and $10^{-16}$, but GENP  
failed on the inputs $\Omega$
both when we used no preprocessing 
and used preprocessing with random 
circulant multipliers of Examples \ref{exc1} and \ref{exunc}. 
In these cases the
  relative residual norms of the output approximations ranged between
$10^{-2}$ and $10^{4}$. In contrast GENP applied to the
inputs preprocessed with Gaussian multipliers 
produced quite reasonable approximations to the solution.
Already after a single
step of iterative refinement,
they have at least matched the level of GEPP. 
Table \ref{tab65a} displays these norms
in some detail.

\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig1.pdf}
	\caption{Norm of $A^{-1}$}	
	\label{fig:fig1}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig2.pdf}
	\caption{Average relative residual norms for GENP by using random multipliers.
	The two broken lines representing one iteration of circulant multipliers are overlapping at the bottom  of the display
	}
	\label{fig:fig2}
\end{figure}

\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig3.pdf}
	\caption{Maximum relative residual norms for GENP by using random multipliers.
	The two broken lines representing one iteration of circulant multipliers are overlapping at the bottom of the display
}
	\label{fig:fig3}
\end{figure}


\subsection{Approximation of the leading singular spaces and 
low-rank appro\-xi\-ma\-tion of 
a matrix}\label{stails}  


We approximated
the $r$-dimensional leading singular spaces of $n\times n$
matrices $A$ that have
numerical rank $r$, and we also
 approximated these matrices with matrices of rank $r$. 
For $n=64,128,256, 512,1024$ and  $r=8,32$ we  generated $n\times n$   random
orthogonal 
matrices $S$ and $T$ and diagonal matrices 
$\Sigma=\diag(\sigma_j)_{j=1}^n$ such that $\sigma_j=1/j,~j=1,\dots,r$,
$\sigma_j=10^{-10},~j=r+1,\dots,n$ (cf. \cite[Section 28.3]{H02}). 
Then we 
computed the input matrices
$A=S_A\Sigma_A T_A^T$, for which $||A||=1$ and $\kappa (A)=10^{10}$.
Furthermore 
we  generated 
 $n\times r$ random matrices $H$ 
and computed the matrices
$B_{r,A}=AH$, $Q_{r,A}=Q(B_{r,A})$, $S_{r,A}$, 
 $T_{r,A}$,  
$Y_{r,A}=Q_{r,A}^TS_{r,A}$, 

 and 
$Q_{r,A}Q_{r,A}^TA$.
Figures \ref{fig:fig4}--\ref{fig:fig7} and
Tables \ref{tab66}--\ref{tab614} display the resulting data on the 
residual norms 
${\rm rn}^{(1)}=||Q_{r,A}Y_{r,A}-S_{r,A}||$ and 
${\rm rn}^{(2)}=||A-Q_{r,A}Q_{r,A}^TA||$, 

obtained in 1000 runs of our tests
for every pair of $n$ and $r$. In these figures and tables 
${\rm rn}^{(1)}$ 
denotes the residual norms 
of the approximations of the matrix bases for the 
leading singular spaces $\mathbb S_{r,A}$, and
${\rm rn}^{(2)}$ denotes the residual norms 
of the approximations
of the matrix $A$ by the rank-$r$ matrix
$Q_{r,A}Q_{r,A}^TA$.

 Figures \ref{fig:fig4} and  \ref{fig:fig5} 
and Tables \ref{tab66}--\ref{tab613} show the norm ${\rm rn}^{(1)}$.  
The last column of each of the tables displays 
the ratio of the observed values ${\rm rn}^{(1)}$ 
 and its upper bound 
$\tilde \Delta_+= \sqrt 2~\frac{\sigma_{r+1}(A)}{\sigma_r(A)}||H||_F ||(T_{r,A}^TH)^{-1}||$ 
 estimated up to the higher order terms (cf. Corollary \ref{coerrb}). 


In our tests we had  $\sigma_{r}(A)=1/r$ and $\sigma_{r+1}(A)=10^{-10}$.
Table \ref{tab66} covers the case where   
we generated  Gaussian multipliers $H$.
Tables \ref{tab610} and \ref{tab613} cover the cases where 
we generated 
random $n\times n$
circulant matrices of Examples \ref{exc1} and \ref{exunc}, respectively,
and applied their $n\times r$ Toeplitz leading blocks as multipliers $H$.

Figures \ref{fig:fig6} and \ref{fig:fig7}  and
Tables \ref{tab67}--\ref{tab614}
show similar results of our tests for the observed residual norms ${\rm rn}^{(2)}$
and their ratios with their upper bounds  $\tilde \Delta'_+=\sigma_{r+1}(A)+2\Delta_+ ||A||$, 
 estimated up to the higher order terms (cf. Corollary \ref{coover}). 

Tables \ref{tab69}--\ref{tab616} show some auxiliary information.
Namely, Table \ref{tab69} displays the data on the ratios $||(T_{r,A}^TH)^{-1}||/||(H_{r,r})^{-1}||$,

where $H_{r,r}$ denotes  the
$r\times r$ leading submatrix of the matrix $H$.


Tables \ref{tab616} and  \ref{tab617} display the average condition numbers 
of Gaussian $n\times n$ matrices and circulant
 $n\times n$  matrices $C$ of Example \ref{exc1},
respectively.
 
The test results are in quite good accordance with our theoretical study
of Gaussian multipliers and
suggest that the power of random circulant and Toeplitz multipliers is  similar 
to the power of Gaussian multipliers,
as in the case of various random structured multipliers
of \cite{HMT11} and \cite{M11}.


\begin{figure}
 \centering
	\includegraphics[width=0.80\textwidth]{fig4.pdf}	
	\caption{Residual norms $rn^{(1)}$ using different random multipliers, case r=8}
	\label{fig:fig4}
\end{figure}


 
\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig5.pdf}
	\caption{Residual norms $rn^{(1)}$ using different random multipliers, case r=32}
	\label{fig:fig5}
\end{figure}


\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig6.pdf}
	\caption{Residual norms  $rn^{(2)}$ using different random multipliers, case r=8}
	\label{fig:fig6}
\end{figure}


 
\begin{figure}
	\centering
		\includegraphics[width=0.80\textwidth]{fig7.pdf}
	\caption{Residual norms $rn^{(2)}$ using different random multipliers, case r=32}
	\label{fig:fig7}
\end{figure}


\section{Conclusions}\label{sconcl}


It is  known that a
 standard Gaussian random matrix
(we refer to it as Gaussian for short) 
has full rank with probability 1 and
is well-con\-di\-tioned with a probability
close to 1.
These properties  motivated our application of 
 Gaussian multipliers
to advancing 
 matrix 
computations.
In particular we preprocessed
 well-con\-di\-tioned nonsingular  input matrices  
by using Gaussian multipliers to support 

GENP (that is, Gaussian elimination with no pivoting) and block
 Gaussian elimination.  
These algorithms can readily fail in 
practical numerical computations without 
preprocessing, but we proved that we can avoid
these problems  
  with a
probability 
close to 1 if we preprocess the input matrix  by pre- or
 post-multiplying it by a
Gaussian matrix.

Our tests were in good accordance with that formal 
result, that is, we generated matrices that 
were hard for GENP, but the problems were consistently 
avoided when we preprocessed the inputs with Gaussian multipliers.
In that case a single loop of iterative refinement was always sufficient
to match the output accuracy of the customary algorithm of  GEPP,
 indicating that GENP with preprocessing 
followed by even a single step of iterative refinement 
is backward stable, similarly to the celebrated result of
\cite{S80}. 

In our tests we observed  similar results even where we applied  Gaussian 
 circulant (rather than Gaussian) multipliers.
Under
this choice we 
 generated only  $n$
random parameters for an $n\times n$
input, and the multiplication stage was
accelerated by a factor of $n/\log (n)$. 
The acceleration factor increases to $n^2/\log (n)$
 when  
the input matrix has the structure of Toeplitz type, but
we could support
numerical stabilization of GENP with Gaussian circulant multipliers  
only empirically.
Moreover, we proved that with a high probability 
Gaussian  circulant multipliers cannot fix numerical instability
 of the elimination algorithms 
for a  specific narrow class of inputs
 (see Theorem \ref{thcgenp}
and Remark  \ref{remext}).

This should motivate the search for 

alternative randomized structured multipliers
 that 
 would be expected to stabilize numerical performance 
of GENP and block Gaussian elimination
for any input or, say, for any 
Toeplitz and Toeplitz-like input matrix.
Among the candidate multipliers,
one can consider the products of random 
circulant and skew-circulant matrices,
possibly used as both pre- and post-multipliers.
Suppose that indeed they are expected to stabilize 
block Gaussian elimination numerically.
Then their support would be valuable for 
numerical application of 
 the MBA celebrated  
algorithm, because it is superfast 
for Toeplitz and Toeplitz-like input matrices
 and hence for their products with 
circulant and skew-circulant matrices 
(cf., e.g.,  \cite{B85},  \cite[Chapter 5]{p01},  and 
\cite{PQZ11}).  
 
We have extended our analysis to the problem of
rank-$r$ approximation of

 an $m\times n$ matrix $A$ 
 having a numerical rank $r$.
With 
 a probability close to $1$
the column set of
the matrix  $AH$,
for  an $n\times (r+p)$ Gaussian matrix $H$
and a small positive
oversampling integer parameter $p$,
approximates a basis for the left  
leading singular
space $ \mathbb S_{r, A}$
associated with the $r$ largest singular values
of an $m\times n$ matrix $A$. 
Having  such an approximate basis 
available, one can readily approximate 
the matrix $A$ by a matrix of rank $r$.

This is an efficient, well developed algorithm
(see \cite{HMT11}), but we 
 proved that this algorithm is expected to
produce a reasonable rank-$r$ approximation 
with Gaussian multipliers already for $p=0$, that is, even 
without customary oversampling,
recommended in \cite{HMT11}.

Then again

in our tests   
the latter techniques 
were  efficient
even where
instead of  Gaussian multipliers
we applied  random Toeplitz 
 multipliers, defined as the maximal leading submatrices of 
random circulant matrices. This has
accelerated the multiplication stage and has
limited
randomization to $n$
 parameters for an $n\times n$
input. 

Formal proof of the power of random structured SRFT 
multipliers with substantial oversampling
is known for low-rank approximation  \cite[Section 11]{HMT11},
and we  immediately extended it to the case when
the products of random unitary circulant 
multipliers and random rectangular permutation matrices 
were applied instead of the SRFT matrices (see Section \ref{srnd}).

A natural research challenge
 is the combination
of our randomized  
multiplicative  preprocessing with
randomized  augmentation and
additive preprocessing, studied 
 in  
\cite{PQ10}, \cite{PQ12},
\cite{PQZC},  \cite{PQZ13}, and \cite{PQZb}.


\bigskip


{\bf Acknowledgements:}
Our research has been supported by NSF Grant CCF--1116736.
 and 
PSC CUNY Awards 64512--0042 and 65792--0043.
We are also grateful to the reviewer 
 for thoughtful helpful comments.
{"config": "arxiv", "file": "1406.5802/000mp_nov_wf.tex"}