remarkRemark assumptionAssumption Randomized Lanczos for null space
computationsDaniel Kressner and Nian Shao

# A randomized small-block Lanczos method for large-scale null space
computations

Daniel Kressner Institute of Mathematics, EPF Lausanne, 1015 Lausanne,
Switzerland<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Nian
Shao11footnotemark: 1

###### Abstract

Computing the null space of a large sparse matrix $A$ is a challenging
computational problem, especially if the nullity—the dimension of the null
space—is large. When using a block Lanczos method for this purpose,
conventional wisdom suggests to use a block size $d$ that is not smaller than
the nullity. In this work, we show how randomness can be utilized to allow for
smaller $d$ without sacrificing convergence or reliability. Even $d=1$,
corresponding to the standard single-vector Lanczos method, becomes a safe
choice. This is achieved by using a small random diagonal perturbation, which
moves the zero eigenvalues of $A^{\mathsf{T}}A$ away from each other, and a
random initial guess. We analyze the effect of the perturbation on the
attainable quality of the null space and derive convergence results that
establish robust convergence for $d=1$. As demonstrated by our numerical
experiments, a smaller block size combined with restarting and partial
reorthogonalization results in reduced memory requirements and computational
effort. It also allows for the incremental computation of the null space,
without requiring a priori knowledge of the nullity.

###### keywords:

Sparse matrix, null space, Krylov subspace method, randomized algorithm

65F15, 65F50, 68W20

##  1 Introduction

This work is concerned with computing the null space of a large sparse matrix
$A\in\mathbb{R}^{m\times n}$, a problem that finds numerous applications. For
example, the connected components of an undirected graph can be identified by
the null space of the graph Laplacian [18]. In computational topology, the
null space of cohomology matrices reflects the generators of the cohomology
groups. The nullity helps determine the Betti numbers of a topological space
[7]. Furthermore, null space calculations play a crucial role in various
engineering problems, such as system dimensionality reduction [34], finite
element models [29], linear equality constrained state estimation [13], grasp
control [25], computer vision [36], and multivariate polynomial root finding
[2].

In several of these applications, the dimension of the null space is not
necessarily small and can easily reach a few hundred or even thousands. For
small- to medium-sized matrices, the singular value decomposition (SVD) [10]
or a suitably modified QR decomposition [8] are the methods of choice for
computing null spaces. For large-sized matrices, these dense methods become
too expensive. There are several structured null space solvers to deal with
specific large-scale situations, such as [3] for banded matrices, [11] for
matrices with small nullity and amenable to sparse LU factorization with
partial pivoting, [9] for sparse matrices amenable to a rank-revealing sparse
QR factorization [4], and [22] for a very tall and skinny matrix.

When $A$ is large and has no particular structure beyond being sparse, one
needs to resort to Krylov subspace methods, which utilize matrix-vector
products with $A$ and its transpose $A^{\mathsf{T}}$. Let $N$ denote the null
space dimension of $A$. When $N=1$, standard choices include the Lanczos
method for computing the smallest eigenvalue and eigenvector of
$A^{\mathsf{T}}A$ or the mathematically equivalent, but numerically more
reliable Golub–Kahan bidiagonalization directly applied to $A$ [10, Chapter
10]. When $N>1$, the Lanczos method (with a single starting vector) becomes
unreliable and prone to miss components of the null space. In fact, in exact
arithmetic only a one-dimensional subspace of the null space can be extracted
via Lanczos [23, Section 12.2]. Block Lanczos methods address this issue by
utilizing a block of $d>1$ starting vectors instead of a single starting
vector [23, Section 13.10]. Common wisdom (see [23, Section 13.12] for
examples) and existing convergence analyses all suggest to choose $d$ not
smaller than the size of the eigenvalue cluster of interest. Applied to null
space computation this would require $d\geq N$ and in exact arithmetic this
condition is indeed necessary. The presence of round-off error blurs the
situation. As an example, consider the $420\times 420$ diagonal matrix

(1) $A=\mathrm{diag}(0,\dotsc,0,1,2,\dotsc,399),$

that is, the first $21$ diagonal entries are zero. When applying block Lanczos
to $A$ in order to compute its null space and choosing a block of $d$
(Gaussian) random starting vectors, one would expect to obtain no more than
$d$ (approximate) zero eigenvalues and corresponding eigenvectors. However,
when executing block Lanczos in double precision arithmetic without breakdown
and declaring eigenvalues less than $10^{-4}$ as zero, we obtained the results
reported in Table 1.

Table 1: Null space dimension obtained when applying block Lanczos with $d$ random starting vectors to the matrix from Eq. 1. Block size $d$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 10 | 12 | 14 | 15 | 20 | 21
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Dimension | 11 | 11 | 16 | 12 | 15 | 15 | 21 | 20 | 21 | 21 | 21 | 21 | 21

While the null space dimension estimated by block Lanczos for small values of
$d$ is erratic and smaller than $21$, it is always larger than $d$. In fact,
already for $d=12$, the correct null space dimension is detected. Experiments
with other matrices lead to similar findings, suggesting that round-off error
breaks some of the curse incurred by eigenvalue clusters, but it does not
fully address the convergence issues of single-vector or small-block Lanczos
methods either.

The situation described above is reminiscent of recent work by Meyer, Musco
and Musco [19] on the single-vector Lanczos method for low-rank approximation.
Taking the randomness of the starting vector into account they establish a
convergence result that still requires the (larger) singular values to be
distinct but the dependence of the complexity on the gaps between eigenvalues
is very mild, in fact logarithmic. Although single-vector or small-block
Krylov subspace methods do not satisfy gap-independent convergence bounds [16,
27, 21, 35], the theory in [19] suggests that small random perturbations of
$A$ can easily break the adverse role of repeated singular values. As we have
seen, perturbations due to round-off are not sufficient to achieve this
effect, but (slightly larger) random perturbations will do, thanks to
eigenvalue repulsion results.

Inspired by [19], we propose a randomized small-block Lanczos method for null
space calculation, outlined in Algorithm 1. Compared to the task of low-rank
approximation, there are two major differences. First, in low-rank
approximation, the largest few singular values are usually unknown and, in
most applications, they do not form tight large clusters. In the context of
null spaces, the desired singular values form one large cluster of zeros.
Second, while convergence to the relevant invariant subspace is not necessary
for low-rank approximation [6], such convergence is imperative to obtain a
good null space approximation.

Input: Matrix $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ and $N$ zero
singular values. Perturbation parameter $\epsilon>0$. Block size $d$. Number
of iterations $\ell$.

Output: Orthonormal basis $V$ approximating null space of $A$.

1 Set $B=A^{\mathsf{T}}A+\epsilon D$, with a diagonal matrix $D$ having
diagonal entries uniformly i.i.d. in $[0,1]$;

2 Draw Gaussian random matrix $\Omega\in\mathbb{R}^{n\times d}$;

3 Perform block Lanczos to compute Lanczos decomposition
$BZ_{\ell}=Z_{\ell}T_{\ell}+Q_{\ell+1}E_{\ell+1}$, where $Z_{\ell}$ is an
orthonormal basis of Krylov subspace
$\mathrm{span}\\{\Omega,B\Omega,\dotsc,B^{\ell-1}\Omega\\}$ and $T_{\ell}$ is
block tridiagonal;

4 Compute an orthonormal basis $V_{Z}$ for the invariant space associated with
$N$ smallest eigenvalues of $T_{\ell}$;

5 return _$V=Z_{\ell}V_{Z}$_ ;

Algorithm 1 Randomized small-block Lanczos for null space computation

######  Remark 1.1.

Line 1 of Algorithm 1 extracts the null space approximation from Ritz vectors
by computing the spectral decomposition of
$T_{\ell}=Z_{\ell}^{\mathsf{T}}BZ_{\ell}$. Alternatively, one could remove the
perturbation and compute the spectral decomposition of
$Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell}$ or, equivalently, the SVD of
$AZ_{\ell}$. Preliminary numerical experiments indicate that this only yields
minor accuracy improvements, at best, not justifying the extra computational
effort needed to perform the alternative procedure.

The rest of this paper is organized as follows.

In Section 2, we analyze the impact of the random perturbation $\epsilon D$
introduced by Algorithm 1. In particular, Theorem 2.2 provides a new
eigenvalue repulsion result for the perturbed zero eigenvalues of
$A^{\mathsf{T}}A$. Our result differs from an existing result [19, Lemma 12]
in two aspects: (1) By focusing on the repulsion of zero eigenvalues rather
than the entire spectrum, we establish $\mathcal{O}(\epsilon)$ gaps instead of
$\mathcal{O}(\epsilon^{2})$ gaps. (2) By considering relative eigenvalue gaps,
we establish lower bound independent of the norm of $A$. At the same time, the
perturbation incurs a limit on the attainable accuracy of the null space
approximation, and we quantify this effect by a perturbation analysis.

Section 3 performs a self-contained convergence analysis of the randomized
single-vector Lanczos method (corresponding to $d=1$ in Algorithm 1). While
our analysis shares similarities with [19], we achieve more refined results by
focusing on $d=1$ and considering the convergence of individual Ritz vectors.
This analysis cumulates in Theorems 3.1 and 3.10, establishing bounds on the
null space approximation, whose sharpness is verified by the numerical
results.

In Section 4, we consider several practical improvements of Algorithm 1,
including preconditioning, restarting, and partial reorthogonalization.

Finally, Section 5 contains several numerical results for applications, such
as the computation of connected graph components and cohomology computation.
These numerical experiments not only validate our theoretical findings
regarding the randomized single-vector Lanczos method but also highlight the
efficiency of the small-block Lanczos method. Notably, it highlights when our
method is preferable compared to other null space solvers, particularly when
the nullity is substantial and memory constraints are a limiting factor.

##### Assumptions and notation

Throughout this paper we assume that $A\in\mathbb{R}^{m\times n}$ satisfies
$m\geq n$ and has $N\geq 1$ zero singular values. In accordance with Step 1 of
Algorithm 1, we define

(2) $B:=A^{\mathsf{T}}A+\epsilon D,$

where $D$ is diagonal with nonnegative diagonal entries and $\epsilon>0$
controls the strength of the perturbation. We let
$\lambda_{1},\ldots,\lambda_{n}$ denote the (nonnegative) eigenvalues of $B$,
sorted such that

$0\leq\lambda_{1}\leq\dotsb\leq\lambda_{n}.$

The initial $\Omega\in\mathbb{R}^{n\times d}$ is chosen as a Gaussian random
matrix, i.e., all entries are i.i.d. (independent and identically distributed)
from the standard normal distribution $\mathcal{N}(0,1)$. We use $\|\cdot\|$
to denote the Euclidean norm of a vector and the spectral norm of a matrix.

######  Remark 1.2.

For null space computation, the cases $m\geq n$ and $m<n$ are fundamentally
different because $m<n$ implies that $A$ has a null space of dimension at
least $n-m$. When $n$ is much larger than $m$, say $n=2m$, one needs to resort
to dense methods like the SVD, in the absence of additional structure.

##  2 Perturbation analysis

The purpose of this section is to quantify the impact of the (random)
perturbation $\epsilon D$ introduced by Algorithm 1.

###  2.1 Perturbation and repulsion of zero eigenvalues

We start by considering the impact of the perturbation on the $N$ zero
eigenvalues of $A^{\mathsf{T}}A$. Suppose that
$0<\epsilon\|D\|<\underline{\sigma}^{2}$, where $\underline{\sigma}$ denotes
the smallest nonzero singular value of $A$. Then the zero and nonzero
eigenvalues of $A^{\mathsf{T}}A$ will not be mixed by the perturbation
$\epsilon D$. Indeed, using that $D$ is positive semi-definite, it follows
from Weyl’s inequality [14, Theorem 4.3.1] that the eigenvalues
$\lambda_{1},\dotsc,\lambda_{n}$ of $B$ satisfy

(3)
$0\leq\lambda_{1}\leq\dotsb\leq\lambda_{N}\leq\epsilon\|D\|<\underline{\sigma}^{2}\leq\lambda_{N+1}\leq\dotsb\leq\lambda_{n}.$

Theorem 2.2 below shows that it is unlikely that the perturbed zero
eigenvalues $\lambda_{1},\dotsc,\lambda_{N}$ remain very close to each other
when $D$ is randomly chosen. To establish this eigenvalue repulsion result, we
follow the proof of [19, Lemma 12] and make use of the following proposition.

######  Proposition 2.1 (Equation (1.11) in [1] attributed to [20]).

Let $B_{0}\in\mathbb{R}^{n\times n}$ be symmetric, and let
$D\in\mathbb{R}^{n\times n}$ be a random diagonal matrix with the diagonal
entries drawn i.i.d. from a distribution with bounded probability density
function $\rho(\cdot)$. Then for any bounded interval
$\mathcal{I}\subset\mathbb{R}$ and $\epsilon>0$, the inequality

$\Pr(B_{0}+\epsilon D\text{ has at least $2$ eigenvalues in }\mathcal{I})\leq
C_{0}(\epsilon^{-1}\lVert\rho\rVert_{\infty}\lvert\mathcal{I}\rvert n)^{2},$

holds for some universal constant $C_{0}$, where $\lvert\mathcal{I}\rvert$ is
the length of $\mathcal{I}$.

######  Theorem 2.2.

Consider the matrix $B$ defined in (2) and suppose that the diagonal entries
of $D$ are drawn i.i.d. from the uniform distribution on $[0,1]$. Let
$\underline{\sigma}$ denote the smallest nonzero singular value of $A$. If
$0<\epsilon<\underline{\sigma}^{2}$, the inequality

$\min_{1\leq i\leq
N-1}\lvert\lambda_{i+1}-\lambda_{i}\rvert\geq\frac{\delta\epsilon}{Cn^{2}},$

holds with probability at least $1-\delta$ for any $0<\delta<1$, where $C$ is
a universal constant.

######  Proof 2.3.

Because of $\|D\|\leq 1$, it follows from (3) that
$\lambda_{1},\dotsc,\lambda_{N}$ are contained in the interval $[0,\epsilon]$.
We cover this interval with overlapping (smaller) intervals
$\mathcal{I}_{j}=[j\gamma,(j+2)\gamma]$ for $j=0,\ldots,\epsilon/\gamma-2$,
where—for the moment—$\gamma$ is an arbitrary number such that
$0<\gamma\leq\epsilon/2$ and $\epsilon/\gamma$ is an integer. If two
neighboring eigenvalues satisfy $\lambda_{i+1}-\lambda_{i}\leq\gamma$, then
there clearly exists an interval $\mathcal{I}_{j}$ such that
$\lambda_{i},\lambda_{i+1}\in\mathcal{I}_{j}$. Thus,

$\displaystyle\Pr\\{\exists i\colon\lambda_{i+1}-\lambda_{i}\leq\gamma\\}$
$\displaystyle\leq$ $\displaystyle\Pr\\{\exists
i,\,j\colon\lambda_{i},\lambda_{i+1}\in\mathcal{I}_{j}\\}\leq\sum_{j=0}^{\epsilon/\gamma-2}\Pr\\{\exists
i\colon\lambda_{i},\lambda_{i+1}\in\mathcal{I}_{j}\\}$ $\displaystyle\leq$
$\displaystyle\sum_{j=0}^{\epsilon/\gamma-2}4C_{0}n^{2}\epsilon^{-2}\gamma^{2}=4C_{0}n^{2}\gamma\epsilon^{-1}\big{(}1-\gamma\epsilon^{-1}\big{)}\leq\frac{4C_{0}n^{2}}{\epsilon/\gamma+1},$

where we used the union bound and Proposition 2.1 in the second and third
inequalities, respectively. The last inequality uses
$x(1-x)=x(1-x^{2})/(1+x)\leq x/(1+x)$ with $x=\gamma\epsilon^{-1}$. The proof
is concluded by setting $C=4C_{0}$,
$\gamma=\epsilon/\lfloor\frac{Cn^{2}}{\delta}\rfloor\geq\frac{\delta\epsilon}{Cn^{2}}$
and using
${\epsilon}/{\gamma}=\lfloor{Cn^{2}}/{\delta}\rfloor>{Cn^{2}}/{\delta}-1$.

###  2.2 Perturbation of null space residual

We will bound the impact of perturbations on the null space residual of $A$
using the following eigenvalue perturbation result.

######  Proposition 2.4.

Consider symmetric matrices $B_{0}$ and $B\in\mathbb{R}^{n\times n}$ with
eigenvalues $0\leq\lambda_{1}^{(0)}\leq\dotsb\leq\lambda^{(0)}_{n}$ and
$0\leq\lambda_{1}\leq\dotsb\leq\lambda_{n}$, respectively. Let $V_{N}$ be an
orthonormal basis for the invariant subspace of $B$ associated with
$\lambda_{1},\dotsc,\lambda_{N}$. Suppose that $B-B_{0}$ is positive semi-
definite, then

$\lVert V_{N}^{\mathsf{T}}B_{0}V_{N}\rVert\leq\lambda_{N}^{(0)}+\lVert
B-B_{0}\rVert.$

If, additionally, $\lVert B-B_{0}\rVert<\lambda_{N+1}^{(0)}-\lambda_{N}^{(0)}$
then

$\lVert V_{N}^{\mathsf{T}}B_{0}V_{N}\rVert\leq\lambda_{N}^{(0)}+\frac{\lVert
B-B_{0}\rVert^{2}}{\lambda_{N+1}^{(0)}-\lambda_{N}^{(0)}-\lVert
B-B_{0}\rVert}.$

######  Proof 2.5.

The first result is a direct consequence of Weyl’s inequality:

$\lVert V_{N}^{\mathsf{T}}B_{0}V_{N}\rVert\leq\lVert
V_{N}^{\mathsf{T}}BV_{N}\rVert=\lambda_{N}\leq\lambda_{N}^{(0)}+\lVert
B-B_{0}\rVert,$

The second result follows from a quadratic residual bound for eigenvalues in
[17, Theorem 1]. To apply this bound, we let $V_{\perp}$ be such that
$[V_{N},V_{\perp}]$ is an $n\times n$ orthogonal matrix and consider the gap

(4) $\displaystyle g$
$\displaystyle=\lambda_{\min}\big{(}V_{\perp}^{\mathsf{T}}B_{0}V_{\perp}\big{)}-\lambda_{N}^{(0)}\geq\lambda_{\min}\big{(}V_{\perp}^{\mathsf{T}}BV_{\perp}\big{)}-\lambda_{N}^{(0)}-\|B-B_{0}\|$
$\displaystyle=\lambda_{N+1}-\lambda_{N}^{(0)}-\|B-B_{0}\|\geq\lambda_{N+1}^{(0)}-\lambda_{N}^{(0)}-\|B-B_{0}\|>0.$

Now, [17, Theorem 1] yields

$\lVert V_{N}^{\mathsf{T}}B_{0}V_{N}\rVert-\lambda_{N}^{(0)}\leq
g^{-1}\lVert(I-V_{N}V_{N}^{\mathsf{T}})B_{0}V_{N}\rVert^{2}\leq g^{-1}\lVert
B-B_{0}\rVert^{2},$

where the second inequality uses that $V_{N}$ spans an invariant subspace of
$B$. Combined with Eq. 4, this proves the desired result.

Proposition 2.4 provides two kinds of perturbation bounds. The first bound
does not require a gap but depends linearly on $\lVert B-B_{0}\rVert$, while
the second one improves this to quadratic dependence at the expense of
requiring a gap between $\lambda_{N}^{(0)}$ and $\lambda_{N+1}^{(0)}$.

######  Theorem 2.6.

Let $V_{N}\in\mathbb{R}^{n\times N}$ be an orthonormal basis of the invariant
subspace associated with the $N$ smallest eigenvalues of
$B=A^{\mathsf{T}}A+\epsilon D$ for a symmetric positive semi-definite matrix
$D$ with $\lVert D\rVert\leq 1$. If $0<\epsilon<\underline{\sigma}^{2}/2$,
where $\underline{\sigma}$ denotes the smallest positive singular value of
$A$, then

$\lVert
AV_{N}\rVert\leq\epsilon(\underline{\sigma}^{2}-\epsilon)^{-1/2}\leq\sqrt{\epsilon}.$

######  Proof 2.7.

The result essentially follows from inserting $B_{0}=A^{\mathsf{T}}A$ into the
second bound of Proposition 2.4.

###  2.3 Principal angles and null space residual

To study the impact of perturbations on subspaces, we recall the definition
and basic properties of principal angles. For later purposes, it will be
helpful to admit subspaces of unequal dimension, as in [10, Section 6.4.3] and
[38].

######  Definition 2.8.

Let the columns of $V_{0}\in\mathbb{R}^{n\times k_{0}}$ and
$V_{1}\in\mathbb{R}^{n\times k_{1}}$ form orthonormal bases of subspaces
$\mathcal{V}_{0}$ and $\mathcal{V}_{1}$, respectively. Let
$k=\min\\{k_{0},k_{1}\\}$ and let
$1\geq\sigma_{1}\geq\sigma_{2}\geq\dotsb\geq\sigma_{k}\geq 0$ denote the
singular values of $V_{0}^{\mathsf{T}}V_{1}$. Then the principal angles
$\angle_{t}(V_{0},V_{1})\in[0,\pi/2]$, $t=1,\dotsc,k$, between
$\mathcal{V}_{0}$ and $\mathcal{V}_{1}$ are defined by

$\cos\angle_{t}(\mathcal{V}_{0},\mathcal{V}_{1})=\sigma_{t}.$

For brevity, we will write $\angle_{t}(\Omega_{0},\Omega_{1})$ instead of
$\angle_{t}(\mathcal{V}_{0},\mathcal{V}_{1})$ for arbitrary, not necessarily
orthonormal, basis for $\mathcal{V}_{0}$ and $\mathcal{V}_{1}$.

Note that the principal angles in Definition 2.8 are sorted in increasing
order, and it is often the largest principal angle $\angle_{k}(V_{0},V_{1})$
that is of interest. In particular, it is well known that
$\sin\angle_{k}(V_{0},V_{1})=\lVert
V_{0}V_{0}^{\mathsf{T}}-V_{1}V_{1}^{\mathsf{T}}\rVert$ holds when
$k_{0}=k_{1}$, implying that the orthogonal projectors are close for small
$\angle_{k}(V_{0},V_{1})$. The other principal angles give partial
information; knowing that some of them are small allows us to conclude that
$\mathcal{V}_{0}$ and $\mathcal{V}_{1}$ nearly share a common subspace. In the
context of null spaces, the following result establishes a residual bound in
terms of principal angles.

######  Proposition 2.9.

Consider the setting of Definition 2.8 and a matrix $A\in\mathbb{R}^{m\times
n}$. Then for every $t=1,\dotsc,k$ there exists a matrix
$W_{t}\in\mathbb{R}^{k_{1}\times t}$ with orthonormal columns such that

$\lVert AV_{1}W_{t}\rVert\leq\lVert AV_{0}\rVert+\lVert
A\rVert\sin\angle_{t}(V_{0},V_{1}).$

######  Proof 2.10.

Let $W_{t}$ be the right singular vectors corresponding to the largest $t$
singular values of $V_{0}^{\mathsf{T}}V_{1}$. Then

$\displaystyle\|(I-V_{0}V_{0}^{\mathsf{T}})V_{1}W_{t}\|^{2}$ $\displaystyle=$
$\displaystyle\max_{\|x\|=1}\|(I-V_{0}V_{0}^{\mathsf{T}})V_{1}W_{t}x\|^{2}$
$\displaystyle=$
$\displaystyle\max_{\|x\|=1}\big{(}\|V_{1}W_{t}x\|^{2}-\|V_{0}V_{0}^{\mathsf{T}}V_{1}W_{t}x\|^{2}\big{)}$
$\displaystyle=$ $\displaystyle
1-\min_{\|x\|=1}\|V_{0}^{\mathsf{T}}V_{1}W_{t}x\|^{2}=1-\sigma_{t}^{2}=\sin^{2}\angle_{t}(V_{0},V_{1}).$

The proof is concluded by the triangular inequality:

$\displaystyle\lVert AV_{1}W_{t}\rVert$ $\displaystyle\leq\lVert
AV_{0}V_{0}^{\mathsf{T}}V_{1}W_{t}\rVert+\lVert
A(I-V_{0}V_{0}^{\mathsf{T}})V_{1}W_{t}\rVert$ $\displaystyle\leq\lVert
AV_{0}\rVert+\lVert A\rVert\sin\angle_{t}(V_{0},V_{1}).$

We conclude this section with a useful characterization of the tangents of
principal angles from [38, Theorem 3.1].

######  Proposition 2.11.

Let $[V_{0},V_{\perp}]$ be an orthonormal matrix, where
$V_{0}\in\mathbb{R}^{n\times k_{0}}$. Suppose
$\Omega_{1}\in\mathbb{R}^{n\times k_{1}}$ has column full rank. If
$\mathrm{rank}(V_{0}^{\mathsf{T}}\Omega_{1})=\mathrm{rank}(\Omega_{1})\leq
k_{0}$, then

$\lVert
V_{\perp}^{\mathsf{T}}\Omega_{1}(V_{0}^{\mathsf{T}}\Omega_{1})^{\dagger}\rVert=\tan\angle_{k_{1}}(V_{0},\Omega_{1}).$

##  3 Convergence analysis

In this section, we study the convergence of Algorithm 1, focusing on the case
of a single random starting vector $\omega$, that is, $d=1$.

###  3.1 Convergence of single-vector Krylov subspaces for small eigenvalues

In the following, we analyze how well the Krylov subspace generated in line 3
of Algorithm 1 contains approximations to (portions of) the invariant subspace
associated with the $N$ smallest eigenvalues of $B$. The analysis mimics
analogous results obtained in [19] for low-rank approximation. Compared to
classical results [28, Chapter VI], our analysis differs by considering the
convergence to a whole subspace and leveraging the randomness of the initial
vector. In the following theorem, the main result of this section;
$\mathsf{Cheb}_{k}(\lambda)$ denotes the Chebyshev polynomial of degree $k$
defined as

(5)
$\mathsf{Cheb}_{k}(\lambda)=\frac{1}{2}\big{(}\big{(}\lambda+\sqrt{\lambda^{2}-1}\big{)}^{k}+\big{(}\lambda+\sqrt{\lambda^{2}-1}\big{)}^{-k}\big{)}.$

######  Theorem 3.1.

For a symmetric matrix $B\in\mathbb{R}^{n\times n}$ with eigenvalues
$\lambda_{1}<\lambda_{2}<\dotsb<\lambda_{N}<\lambda_{N+1}\leq\dotsb\leq\lambda_{n}$,
let $V_{N}\in\mathbb{R}^{n\times N}$ be an orthonormal basis of the invariant
subspace associated with $\lambda_{1},\dotsc,\lambda_{N}$. Choose integers
$t,\ell$ such that $1\leq t\leq N$ and $t\leq\ell\leq n$. Let
$K_{\ell}=[\omega,B\omega,\dotsc,B^{\ell-1}\omega]$, where
$\omega\in\mathbb{R}^{n}$ is a Gaussian random vector. Then, for any
$0<\delta<1$, the inequality

(6)
$\tan\angle_{t}(V_{N},K_{\ell})\leq\frac{N\sqrt{5\pi(n-N)\log(n/\delta)}}{\sqrt{2}\delta(t-1)!\cdot\mathsf{Cheb}_{\ell-t}\bigl{(}1+2\frac{\lambda_{N+1}-\lambda_{N}}{\lambda_{n}-\lambda_{N+1}}\bigr{)}}\Bigl{(}\frac{2}{\mathsf{relgap}}\Bigr{)}^{t-1}$

holds with probability at least $1-2\delta$, where

(7) $\mathsf{relgap}:=\min_{1\leq i\leq
N-1}(\lambda_{i+1}-\lambda_{i})/(\lambda_{n}-\lambda_{1}).$

Combining Theorem 3.1 with Theorem 2.2, we obtain the following result.

######  Corollary 3.2.

Suppose that the matrix $B$ in Theorem 3.1 takes the form
$B=A^{\mathsf{T}}A+\epsilon D$, where $A$ is an $m\times n$ matrix with null
space dimension $N$ and $D$ is a random diagonal matrix with the diagonal
entries drawn i.i.d. from the uniform distribution on $[0,1]$. Let
$\underline{\sigma}$ and $\overline{\sigma}$ be the smallest and largest
nonzero singular values of $A$, respectively. If
$0<\epsilon<\underline{\sigma}^{2}$, then

(8)
$\tan\angle_{t}(V_{N},K_{\ell})\leq\frac{N\sqrt{(n-N)\log(n/\delta)}\bigl{(}C(n^{2}\overline{\sigma}^{2})/(\delta\epsilon)\bigr{)}^{t-1}}{\delta\cdot(t-1)!\cdot\mathsf{Cheb}_{\ell-t}\bigl{(}1+\frac{2(\underline{\sigma}^{2}-\epsilon)}{\overline{\sigma}^{2}-\underline{\sigma}^{2}+\epsilon}\bigr{)}}$

holds with probability at least $1-3\delta$, where $C$ is a universal
constant.

The rest of this section is concerned with the proof of Theorem 3.1, which
follows the ideas from [19] and is divided into two parts. First, in Section
3.1.1, we will establish Eq. 6 for $\ell=t$, that is,

(9)
$\tan\angle_{t}(V_{N},K_{t})\leq\frac{N\sqrt{5\pi(n-N)\log(n/\delta)}}{\sqrt{2}\delta(t-1)!}\Bigl{(}\frac{2}{\mathsf{relgap}}\Bigr{)}^{t-1}.$

Second, in Section 3.1.2, we prove

(10)
$\tan\angle_{t}(V_{N},K_{\ell})\leq\frac{\tan\angle_{t}(V_{N},K_{t})}{\mathsf{Cheb}_{\ell-t}\big{(}1+2\frac{\lambda_{N+1}-\lambda_{N}}{\lambda_{n}-\lambda_{N+1}}\big{)}}.$

Plugging Eq. 9 into Eq. 10 then proves the inequality Eq. 6 of Theorem 3.1.

####  3.1.1 The case $\ell=t$. Proof of Eq. 9

Following [19, Theorem 3], the ability of a (single-vector) Krylov subspace to
capture a subspace is closely related to properties of Vandermonde matrices.

######  Lemma 3.3.

For $\lambda_{1}<\dotsb<\lambda_{N}\leq\lambda_{N+1}\leq\dotsb\leq\lambda_{n}$
and $1\leq t\leq N$, consider the Vandermonde matrices

(11)
$\mathsf{Van}_{N}=\begin{bmatrix}1&\lambda_{1}&\ldots&\lambda_{1}^{t-1}\\\
1&\lambda_{2}&\ldots&\lambda_{2}^{t-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
1&\lambda_{N}&\ldots&\lambda_{N}^{t-1}\end{bmatrix},\quad\mathsf{Van}_{\perp}=\begin{bmatrix}1&\lambda_{N+1}&\ldots&\lambda_{N+1}^{t-1}\\\
1&\lambda_{N+2}&\ldots&\lambda_{N+2}^{t-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
1&\lambda_{n}&\ldots&\lambda_{n}^{t-1}\end{bmatrix}.$

Then

$\lVert\mathsf{Van}_{\perp}\mathsf{Van}_{N}^{\dagger}\rVert\leq\frac{\sqrt{n-N}}{(t-1)!}\Bigl{(}\frac{2}{\mathsf{relgap}}\Bigr{)}^{t-1}$

with $\mathsf{relgap}$ defined as in Eq. 7.

######  Proof 3.4.

For an arbitrary vector
$x=[x_{1},\dotsc,x_{t}]^{\mathsf{T}}\in\mathbb{R}^{t}$, we have

$\lVert\mathsf{Van}_{N}x\rVert^{2}=\sum_{i=1}^{N}p^{2}(\lambda_{i})\quad\text{and}\quad\lVert\mathsf{Van}_{\perp}x\rVert^{2}=\sum_{j=N+1}^{n}p^{2}(\lambda_{j}),$

with the polynomial $p(\lambda)=\sum_{k=1}^{t}x_{k}\lambda^{k-1}$. By the
assumptions, $\mathsf{Van}_{N}$ has full column rank and

(12)
$\lVert\mathsf{Van}_{\perp}\mathsf{Van}_{N}^{\dagger}\rVert^{2}=\max_{x\in\mathbb{R}^{t}}\frac{\lVert\mathsf{Van}_{\perp}x\rVert^{2}}{\lVert\mathsf{Van}_{N}x\rVert^{2}}=\max_{p\in\mathcal{P}_{t-1}}\frac{\sum_{j=N+1}^{n}p^{2}(\lambda_{j})}{\sum_{i=1}^{N}p^{2}(\lambda_{i})},$

where $\mathcal{P}_{t-1}$ denotes the vector space containing all polynomials
of degree at most $t-1$. Choosing the Lagrange basis

$p_{i}(\lambda)=\prod_{\begin{subarray}{c}1\leq k\leq t\\\ k\neq
i\end{subarray}}\frac{\lambda-\lambda_{k}}{\lambda_{i}-\lambda_{k}}\in\mathcal{P}_{t-1},\quad
i=1,\dotsc,t,$

we can express an arbitrary polynomial $p\in\mathcal{P}_{t-1}$ as
$p(\lambda)=p(\lambda_{1})p_{1}(\lambda)+\dotsb+p(\lambda_{t})p_{t}(\lambda).$
By the Cauchy–Schwarz inequality,

(13)
$p^{2}(\lambda)\leq\sum_{i=1}^{t}p^{2}(\lambda_{i})\sum_{i=1}^{t}p_{i}^{2}(\lambda)\leq\sum_{i=1}^{N}p^{2}(\lambda_{i})\Bigl{(}\sum_{i=1}^{t}\lvert
p_{i}(\lambda)\rvert\Bigr{)}^{2}.$

To process the last term, we assume that
$\lambda_{N+1}\leq\lambda\leq\lambda_{n}$, which implies
$\lvert\lambda_{i}-\lambda_{k}\rvert/(\lambda-\lambda_{k})\geq\lvert
i-k\rvert\cdot\mathsf{relgap}$ for $1\leq i\neq k\leq t$. In turn,

$\sum_{i=1}^{t}\lvert
p_{i}(\lambda)\rvert=\sum_{i=1}^{t}\prod_{\begin{subarray}{c}1\leq k\leq t\\\
k\neq
i\end{subarray}}\frac{\lambda-\lambda_{k}}{\lvert\lambda_{i}-\lambda_{k}\rvert}\leq\frac{1}{\mathsf{relgap}^{t-1}}\sum_{i=1}^{t}\prod_{\begin{subarray}{c}1\leq
k\leq t\\\ k\neq i\end{subarray}}\frac{1}{\lvert
i-k\rvert}=\Bigl{(}\frac{2}{\mathsf{relgap}}\Bigr{)}^{t-1}\frac{1}{(t-1)!},$

where the last equality utilizes the combinatorial identity

$\sum_{i=1}^{t}\prod_{\begin{subarray}{c}1\leq k\leq t\\\ k\neq
i\end{subarray}}\frac{1}{\lvert
i-k\rvert}=\sum_{i=1}^{t}\frac{1}{(i-1)!(t-i)!}=\frac{1}{(t-1)!}\sum_{i=1}^{t}\binom{t-1}{i-1}=\frac{2^{t-1}}{(t-1)!}.$

Plugging this bound into (13) for $\lambda=\lambda_{j}$ gives

$\displaystyle\frac{\sum_{j=N+1}^{n}p^{2}(\lambda_{j})}{\sum_{i=1}^{N}p^{2}(\lambda_{i})}\leq\frac{n-N}{\bigl{(}(t-1)!\bigr{)}^{2}}\Bigl{(}\frac{2}{\mathsf{relgap}}\Bigr{)}^{2(t-1)}.$

Combined with Eq. 12, this completes the proof.

By leveraging the randomness of $\omega$, Proposition 2.11 allows us to prove
Eq. 9.

######  Proof 3.5 (Proof of Eq. 9).

Because $V_{N}^{\mathsf{T}}\omega$ is a Gaussian random vector, the matrix
$\mathrm{diag}(V_{N}^{\mathsf{T}}\omega)$ is invertible (almost surely). Using
$V_{N}^{\mathsf{T}}K_{t}=\mathrm{diag}(V_{N}^{\mathsf{T}}\omega)\mathsf{Van}_{N}$,
the invertibility of this matrix allows us to conclude that the Krylov basis
$K_{t}=[\omega,B\omega,\dotsc,B^{t-1}\omega]$ has full rank:

(14)
$\mathrm{rank}(K_{t})\geq\mathrm{rank}(V_{N}^{\mathsf{T}}K_{t})=\mathrm{rank}(\mathsf{Van}_{N})=t,$

where the last equality utilizes that the $N$ smallest eigenvalues of $B$ are
mutually different. Letting $V_{\perp}$ denote an orthonormal basis of the
invariant subspace associated with $\lambda_{N+1},\dotsc,\lambda_{n}$, we now
apply Proposition 2.11 to estimate

(15) $\displaystyle\tan^{2}\angle_{t}(V_{N},K_{t})$ $\displaystyle=\lVert
V_{\perp}^{\mathsf{T}}K_{t}(V_{N}^{\mathsf{T}}K_{t})^{\dagger}\rVert^{2}$
$\displaystyle=\lVert\mathrm{diag}(V_{\perp}^{\mathsf{T}}\omega)\mathsf{Van}_{\perp}(\mathrm{diag}(V_{N}^{\mathsf{T}}\omega)\mathsf{Van}_{N})^{\dagger}\rVert^{2}$
$\displaystyle\leq\frac{\sigma_{\max}^{2}(\mathrm{diag}(V_{\perp}^{\mathsf{T}}\omega))}{\sigma_{\min}^{2}(\mathrm{diag}(V_{N}^{\mathsf{T}}\omega))}\lVert\mathsf{Van}_{\perp}\mathsf{Van}_{N}^{\dagger}\rVert^{2}.$

Because $V_{\perp}^{\mathsf{T}}\omega$ and $V_{N}^{\mathsf{T}}\omega$ are
independent Gaussian random vectors, we can use concentration and anti-
concentration results for chi-squared random variables from [15, Lemma 1] and
[19, Lemma 4], respectively, to conclude that

$\frac{\sigma_{\max}^{2}(\mathrm{diag}(V_{\perp}^{\mathsf{T}}\omega))}{\sigma_{\min}^{2}(\mathrm{diag}(V_{N}^{\mathsf{T}}\omega))}\leq\frac{5\pi
N^{2}\log\bigl{(}n/\delta\bigr{)}}{2\delta^{2}}$

holds with probability at least $1-2\delta$. Combined with the bound on
$\lVert\mathsf{Van}_{\perp}\mathsf{Van}_{N}^{\dagger}\rVert$ from Lemma 3.3,
this completes the proof of Eq. 9.

####  3.1.2 The case $\ell>t$. Proof of Eq. 10

Our proof crucially relies on a simple but important observation from [19,
(1.3)], which establishes the following equivalence between the single-vector
Krylov subspace $\mathrm{span}\\{K_{\ell}\\}$ and a block Krylov subspace
generated from a specific starting matrix:

(16)
$\mathrm{span}\\{K_{\ell}\\}=\mathrm{span}\\{\omega,B\omega,\dotsc,B^{\ell-1}\omega\\}=\mathrm{span}\\{K_{t},BK_{t},\dotsc,B^{\ell-t}K_{t}\\},$

with $t\leq\ell\leq n$. The following lemma is a straightforward extension of
known convergence results for Krylov subspaces to block Krylov subspaces.

######  Lemma 3.6.

For a symmetric matrix $B\in\mathbb{R}^{n\times n}$ with eigenvalues
$\lambda_{1}<\lambda_{2}<\dotsb<\lambda_{N}<\lambda_{N+1}\leq\dotsb\leq\lambda_{n}$,
let $V_{N}\in\mathbb{R}^{n\times N}$ be an orthonormal basis of the invariant
subspace associated with $\lambda_{1},\dotsc,\lambda_{N}$. Choose integers
$t,\ell$ such that $1\leq t\leq N$ and $t\leq\ell\leq n$. For an $n\times t$
matrix $\Omega$ such that $\mathrm{rank}(V_{N}^{\mathsf{T}}\Omega)=t$, let
$K=[\Omega,B\Omega,\dotsc,B^{\ell-t}\Omega]\in\mathbb{R}^{n\times(\ell-t+1)t}$.
Then

(17)
$\tan\angle_{t}(V_{N},K)\leq\frac{\tan\angle_{t}(V_{N},\Omega)}{\mathsf{Cheb}_{\ell-t}\Bigl{(}1+2\frac{\lambda_{N+1}-\lambda_{N}}{\lambda_{n}-\lambda_{N+1}}\Bigr{)}},$

where $\mathsf{Cheb}_{\ell-t}$ is the Chebyshev polynomial Eq. 5 of degree
$\ell-t$.

######  Proof 3.7.

Let
$p(\lambda)=\alpha_{0}+\alpha_{1}\lambda+\dotsb\alpha_{\ell-t}\lambda^{\ell-t}$
be such that $p(\lambda)\not=0$ for $\lambda_{1}\leq\lambda\leq\lambda_{N}$.
Defining the matrix
$H_{K}=[\alpha_{0}I_{t},\alpha_{1}I_{t},\dotsc,\alpha_{\ell-t}I_{t}]^{\mathsf{T}}$,
we have that

$V_{N}^{\mathsf{T}}KH_{K}=\sum_{i=0}^{\ell-t}\alpha_{i}V_{N}^{\mathsf{T}}B^{i}\Omega=\sum_{i=0}^{\ell-t}\alpha_{i}\Lambda_{N}^{i}V_{N}^{\mathsf{T}}\Omega=p(\Lambda_{N})V_{N}^{\mathsf{T}}\Omega,$

and an analogous formula holds when $V_{N}$ is replaced by $V_{\perp}$. By
assumption, $p(\Lambda_{N})$ is invertible and thus $V_{N}^{\mathsf{T}}KH_{K}$
has rank $t$. Additionally, using that
$\mathrm{span}\\{KH_{K}\\}\subset\mathrm{span}\\{K\\}$ and the expression for
$\tan\angle_{t}$ from Proposition 2.11, we obtain that

$\displaystyle\tan\angle_{t}(V_{N},K)$
$\displaystyle\leq\tan\angle_{t}(V_{N},KH_{K})=\lVert
V_{\perp}^{\mathsf{T}}KH_{K}(V_{N}^{\mathsf{T}}KH_{K})^{\dagger}\rVert$
$\displaystyle=\lVert
p(\Lambda_{\perp})V_{\perp}^{\mathsf{T}}\Omega(V_{N}^{\mathsf{T}}\Omega)^{\dagger}p(\Lambda_{N})^{-1}\rVert$
$\displaystyle\leq\frac{\max\limits_{N+1\leq j\leq n}\lvert
p(\lambda_{j})\rvert}{\min\limits_{1\leq i\leq N}\lvert
p(\lambda_{i})\rvert}\lVert
V_{\perp}^{\mathsf{T}}\Omega(V_{N}^{\mathsf{T}}\Omega)^{\dagger}\rVert=\frac{\max\limits_{N+1\leq
j\leq n}\lvert p(\lambda_{j})\rvert}{\min\limits_{1\leq i\leq N}\lvert
p(\lambda_{i})\rvert}\tan\angle_{t}(V_{N},\Omega).$

The proof is completed by setting
$p(\lambda)=\mathsf{Cheb}_{\ell-t}\bigl{(}1+\frac{2(\lambda_{N+1}-\lambda)}{\lambda_{n}-\lambda_{N+1}}\bigr{)}$
and using that a Chebyshev polynomial is bounded by one in magnitude inside
$[-1,1]$ and monotone outside that interval.

###  3.2 Convergence of Algorithm 1

In this subsection, we investigate the quality of the approximate null space
produced by Algorithm 1 via applying the Rayleigh–Ritz procedure to the Krylov
subspace. Unlike existing results on the block Lanczos method [16], we also
present the convergence of lower-dimensional subspaces by bounding the
principal angles $\angle_{t}$ for all $1\leq t\leq N$. We start with a result
that establishes the quasi-optimality of Rayleigh–Ritz.

######  Lemma 3.8.

Under the setting of Theorem 2.6, let $V$ be the approximate null space basis
returned after $\ell\geq N$ iterations of Algorithm 1. Let
$\tilde{\lambda}_{t}^{(0)}$ be the $t$th smallest eigenvalue of
$Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell}$, where $1\leq t\leq\ell$ and
$Z_{\ell}$ is an orthonormal basis of
$\mathrm{span}\\{\omega,B\omega,\dotsc,B^{\ell-1}\omega\\}$. Then for all
$1\leq t\leq N$, there exists a matrix $W_{t}\in\mathbb{R}^{N\times t}$ with
orthonormal columns such that

(18) $\lVert AVW_{t}\rVert^{2}\leq\tilde{\lambda}_{t}^{(0)}+\epsilon.$

Furthermore, if $t=N$ and
$0<\epsilon<\underline{\sigma}^{2}-\tilde{\lambda}_{N}^{(0)}$, then

(19) $\lVert
AV\rVert^{2}\leq\tilde{\lambda}_{N}^{(0)}+\frac{\epsilon^{2}}{\underline{\sigma}^{2}-\tilde{\lambda}_{N}^{(0)}-\epsilon}.$

######  Proof 3.9.

Let $V_{Z}$ and $V_{Z}W_{t}$ be two matrices with orthonormal columns, where
their columns form orthonormal bases for the invariant subspaces associated
with the $N$ and $t$ smallest eigenvalues of $Z_{\ell}^{\mathsf{T}}BZ_{\ell}$,
respectively. By definition, $V=Z_{\ell}V_{Z}$. The gap-independent
perturbation result in Proposition 2.4 allows us to conclude (18):

$\lVert
AVW_{t}\rVert^{2}=\lVert(V_{Z}W_{t})^{\mathsf{T}}(Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell})(V_{Z}W_{t})\rVert\leq\tilde{\lambda}_{t}^{(0)}+\lVert
Z_{\ell}^{\mathsf{T}}(B-A^{\mathsf{T}}A)Z_{\ell}\rVert\leq\tilde{\lambda}_{t}^{(0)}+\epsilon,$

where we used that $\tilde{\lambda}_{t}^{(0)}$ is the $t$th smallest
eigenvalue of $Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell}$.

For Eq. 19, the result when $\ell=N$ is straightforward since $\lVert
AV\rVert^{2}=\tilde{\lambda}_{N}^{(0)}$. To show the case $\ell>N$, we first
note that eigenvalue interlacing implies
$\tilde{\lambda}_{N+1}^{(0)}\geq\underline{\sigma}^{2}$. Because of
$\|Z_{\ell}^{\mathsf{T}}BZ_{\ell}-Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell}\|\leq\epsilon<\underline{\sigma}^{2}-\tilde{\lambda}_{N}^{(0)}\leq\tilde{\lambda}_{N+1}^{(0)}-\tilde{\lambda}_{N}^{(0)}$,
we can apply the second, gap-dependent result in Proposition 2.4 to establish
(19).

By combining the eigenvalue repulsion result from Theorem 2.2, the Krylov
subspace convergence result from Theorem 3.1 and the quasi-optimality of
Rayleigh–Ritz from Lemma 3.8, we can conclude our main convergence result on
the convergence of Algorithm 1.

######  Theorem 3.10.

Let $\overline{\sigma}$ and $\underline{\sigma}$ denote the largest and
smallest nonzero singular values of a matrix $A\in\mathbb{R}^{m\times n}$
having null space dimension $N$. Let $V$ be the approximate null space basis
returned by Algorithm 1 with $0<\epsilon<\underline{\sigma}^{2}/3$ and
$N\leq\ell\leq n$. For every $1\leq t\leq N$, there exists a matrix
$W_{t}\in\mathbb{R}^{N\times t}$ with orthonormal columns such that

(20) $\lVert
AVW_{t}\rVert^{2}\leq\Bigl{(}\frac{\epsilon}{(\underline{\sigma}^{2}-\epsilon)^{1/2}}+\overline{\sigma}\tan\angle_{t}(V_{N},Z_{\ell})\Bigr{)}^{2}+\epsilon,$

holds with probability at least $1-3\delta$, where
$\tan\angle_{t}(V_{N},Z_{\ell})$ is bounded according to Eq. 8. For $t=N$,
when

(21)
$\tan\angle_{N}(V_{N},Z_{\ell})\leq\frac{2\underline{\sigma}^{2}-5\epsilon}{3\overline{\sigma}(\underline{\sigma}^{2}-\epsilon)^{1/2}},$

the bound can be improved to

(22) $\lVert
AV\rVert^{2}\leq\Bigl{(}\frac{\epsilon}{(\underline{\sigma}^{2}-\epsilon)^{1/2}}+\overline{\sigma}\tan\angle_{N}(V_{N},Z_{\ell})\Bigr{)}^{2}+\frac{9\epsilon^{2}}{5(\underline{\sigma}^{2}-\epsilon)}.$

######  Proof 3.11.

Using Lemma 3.8, it suffices to bound $\tilde{\lambda}_{t}^{(0)}$, the $t$th
smallest eigenvalue of $Z_{\ell}^{\mathsf{T}}A^{\mathsf{T}}AZ_{\ell}$. Using
Proposition 2.9, Theorem 2.6, and
$\lvert\sin\angle\rvert\leq\lvert\tan\angle\rvert$, we derive

$\sqrt{\tilde{\lambda}_{t}^{(0)}}\leq\lVert AV_{N}\rVert+\lVert
A\rVert\sin\angle_{t}(V_{N},Z_{\ell})\leq\frac{\epsilon}{(\underline{\sigma}^{2}-\epsilon)^{1/2}}+\overline{\sigma}\tan\angle_{t}(V_{N},Z_{\ell}).$

Combined with Eq. 18, this gives Eq. 20.

In order to prove Eq. 22, we first note that the condition Eq. 21 implies

$\underline{\sigma}^{2}-\epsilon-\tilde{\lambda}_{N}^{(0)}\geq\underline{\sigma}^{2}-\epsilon-\Bigl{(}\frac{\epsilon}{(\underline{\sigma}^{2}-\epsilon)^{1/2}}+\overline{\sigma}\tan\angle_{N}(V_{N},Z_{\ell})\Bigr{)}^{2}\geq\frac{5(\underline{\sigma}^{2}-\epsilon)}{9}>0.$

Thus, the condition for the bound Eq. 19 in Lemma 3.8 is met, which implies
the desired result Eq. 22.

######  Remark 3.12.

Let $\kappa=\overline{\sigma}/\underline{\sigma}$ be the effective condition
number of $A$, and $\epsilon=c\underline{\sigma}^{2}$, where $0<c<1/3$.
Theorem 3.10 states that, with probability at least $1-3\delta$,

$\ell=\mathcal{O}\biggl{(}1+\frac{\log\bigl{(}(n\kappa)/(t\delta\epsilon)\bigr{)}}{\log(1+2\kappa^{-1})}\biggr{)}t$

iterations of Algorithm 1 are needed to partially approximate the null space
of $A$ in the sense of

$\tan\angle_{t}(V_{N},Z_{\ell})\leq\frac{\epsilon}{\sqrt{6}\overline{\sigma}\underline{\sigma}}\leq\frac{2\underline{\sigma}^{2}-5\epsilon}{3\overline{\sigma}(\underline{\sigma}^{2}-\epsilon)^{1/2}},\quad\text{and}\quad\lVert
AVW_{t}\rVert^{2}\leq\frac{17\epsilon}{9},$

where we used Stirling’s formula and a known bound for Chebyshev polynomials
in [28, Section 4.4.1]. Let $BZ_{\ell}=Z_{\ell}T_{\ell}+Q_{\ell+1}E_{\ell+1}$
represent the Lanczos decomposition computed by Algorithm 1. Then,
$T_{\ell}-3\epsilon I$ has at least $t$ negative eigenvalues, all originating
from Ritz vectors close to the null space. This assertion stems from the fact
that the $(N+1)$st smallest eigenvalue of $B$ is at least
$\underline{\sigma}^{2}=\epsilon/c>3\epsilon$. Moreover, for $t=N$, the error
bound improves to

$\lVert
AV\rVert^{2}\leq\Bigl{(}\frac{\epsilon}{(\underline{\sigma}^{2}-\epsilon)^{1/2}}+\frac{\epsilon}{\sqrt{6}\underline{\sigma}}\Bigr{)}^{2}+\frac{9\epsilon^{2}}{5(\underline{\sigma}^{2}-\epsilon)}\leq\frac{161\epsilon^{2}}{30\underline{\sigma}^{2}}.$

Our theoretical results highlight a crucial distinction between the
convergence behaviors of the single-vector and block Lanczos methods. In
contrast to the block Lanczos method [16], where Ritz vectors converge
collectively, the single-vector Lanczos method exhibits individual convergence
of Ritz vectors. This individual convergence allows for deflation and adaptive
detection of the null space dimension, which is particularly advantageous in
the typical scenario where prior information about $N$ is unavailable.

##  4 Practical techniques for implementation

In this section, we discuss several improvements and implementation details to
make Algorithm 1 competitive.

###  4.1 Restarting

Algorithm 1 needs to store the orthonormal basis $Z_{\ell}$ of the Krylov
subspace to allow for reorthogonalization (see below) and Ritz vector
computations. As the number of iterations $\ell$ required for convergence
increases, storing $Z_{\ell}$ becomes challenging due to memory limitations.
Additionally, the cost of Rayleigh–Ritz extraction grows significantly with
$\ell$. To address these issues, we employ a restart mechanism, akin to the
Krylov–Schur method [32, 37].

To conveniently describe our restarting mechanism, we use the superscripts
$\cdot^{(\mathrm{old})}$ and $\cdot^{(\mathrm{new})}$ to denote variables
before and after restarting, respectively. We assume that the dimension of the
Krylov subspace before restarting is $\ell d$ and that $\ell_{0}d$ Ritz pairs
are retained after restarting, where $\ell_{0}<\ell$ and $d$ represents the
block size of the initial matrix $\Omega$.

Suppose we have a Lanczos decomposition of the form

(23)
$BZ_{\ell}^{(\mathrm{old})}=Z_{\ell}^{(\mathrm{old})}T_{\ell}^{(\mathrm{old})}+Q_{\ell+1}^{(\mathrm{old})}E_{\ell+1}^{(\mathrm{old})},$

where $[Z_{\ell}^{(\mathrm{old})},Q_{\ell+1}^{(\mathrm{old})}]$ is a matrix
with orthonormal columns,
$E_{\ell+1}^{(\mathrm{old})}=[0,\Gamma_{\ell}^{(\mathrm{old})}]$, and

(24)
$T_{\ell}^{(\mathrm{old})}=\begin{bmatrix}\Theta^{(\mathrm{old})}&F^{(\mathrm{old})}&&&&\\\
(F^{(\mathrm{old})})^{\mathsf{T}}&\Xi_{\ell_{0}+1}^{(\mathrm{old})}&\Gamma_{\ell_{0}+1}^{(\mathrm{old})}&&&\\\
&(\Gamma_{\ell_{0}+1}^{(\mathrm{old})})^{\mathsf{T}}&\Xi_{\ell_{0}+2}^{(\mathrm{old})}&\Gamma_{\ell_{0}+2}^{(\mathrm{old})}&&\\\
&&\ddots&\ddots&\ddots&\\\
&&&(\Gamma_{\ell-2}^{(\mathrm{old})})^{\mathsf{T}}&\Xi_{\ell-1}^{(\mathrm{old})}&\Gamma_{\ell-1}^{(\mathrm{old})}\\\
&&&&(\Gamma_{\ell-1}^{(\mathrm{old})})^{\mathsf{T}}&\Xi_{\ell}^{(\mathrm{old})}\\\
\end{bmatrix}$

is a block tridiagonal matrix, partitioned such that
$\Theta^{(\mathrm{old})}\in\mathbb{R}^{\ell_{0}d\times\ell_{0}d}$ contains the
information from the previous cycle. To restart and keep $\ell_{0}d$ Ritz
vectors, we apply the Rayleigh–Ritz procedure:

$T_{\ell}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})}=V_{Z}^{(\mathrm{new})}\Theta^{(\mathrm{new})}\quad\text{and}\quad
Z_{\ell_{0}}^{(\mathrm{new})}=Z_{\ell}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})}\in\mathbb{R}^{n\times\ell_{0}d},$

where $V_{Z}^{(\mathrm{new})}\in\mathbb{R}^{\ell d\times\ell_{0}d}$ is an
orthonormal basis of the invariant subspace corresponding to the $\ell_{0}d$
smallest eigenvalues of $T_{\ell}^{(\mathrm{old})}$. Multiplying
$V_{Z}^{(\mathrm{new})}$ from the right to the Lanczos decomposition Eq. 23
leads to

$BZ_{\ell_{0}}^{(\mathrm{new})}=Z_{\ell_{0}}^{(\mathrm{new})}\Theta^{(\mathrm{new})}+Q_{\ell+1}^{(\mathrm{old})}E_{\ell+1}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})}.$

We restart the block Lanczos method by defining the quantities

$\displaystyle
Q_{\ell_{0}+1}^{(\mathrm{new})}=Q_{\ell+1}^{(\mathrm{old})}\quad\text{and}\quad
Z_{\ell_{0}+1}^{(\mathrm{new})}=[Z_{\ell_{0}}^{(\mathrm{new})},Q_{\ell_{0}+1}^{(\mathrm{new})}]\in\mathbb{R}^{n\times(\ell_{0}+1)d},$
$\displaystyle
F^{(\mathrm{new})}=(E_{\ell+1}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})})^{\mathsf{T}}\quad\text{and}\quad\Xi_{\ell_{0}+1}^{(\mathrm{new})}=(Q_{\ell_{0}+1}^{(\mathrm{new})})^{\mathsf{T}}BQ_{\ell_{0}+1}^{(\mathrm{new})},$
$\displaystyle\bigl{(}I-Z_{\ell_{0}+1}^{(\mathrm{new})}(Z_{\ell_{0}+1}^{(\mathrm{new})})^{\mathsf{T}}\bigr{)}BQ_{\ell_{0}+1}^{(\mathrm{new})}=Q_{\ell_{0}+2}^{(\mathrm{new})}\Gamma_{\ell_{0}+1}^{(\mathrm{new})},$

where the last line denotes a QR decomposition. After restarting, the new
Lanczos decomposition takes the form

$\displaystyle BZ_{\ell_{0}+1}^{(\mathrm{new})}$
$\displaystyle=[BZ_{\ell_{0}}^{(\mathrm{new})},BQ_{\ell_{0}+1}^{(\mathrm{new})}]$
$\displaystyle=[Z_{\ell_{0}}^{(\mathrm{new})},Q_{\ell_{0}+1}^{(\mathrm{new})}]\begin{bmatrix}\Theta^{(\mathrm{new})}&(E_{\ell+1}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})})^{\mathsf{T}}\\\
E_{\ell+1}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})}&\Xi_{\ell_{0}+1}^{(\mathrm{new})}\end{bmatrix}+Q_{\ell_{0}+2}^{(\mathrm{new})}[O,\Gamma_{\ell_{0}+1}^{(\mathrm{new})}]$
$\displaystyle=Z_{\ell_{0}+1}^{(\mathrm{new})}\begin{bmatrix}\Theta^{(\mathrm{new})}&F^{(\mathrm{new})}\\\
(F^{(\mathrm{new})})^{\mathsf{T}}&\Xi_{\ell_{0}+1}^{(\mathrm{new})}\end{bmatrix}+Q_{\ell_{0}+2}^{(\mathrm{new})}E_{\ell_{0}+2}^{(\mathrm{new})},$

where
$E_{\ell_{0}+2}^{(\mathrm{new})}=[O,\Gamma_{\ell_{0}+1}^{(\mathrm{new})}]$,
and the relationship

$E_{\ell+1}^{(\mathrm{old})}V_{Z}^{(\mathrm{new})}=(Q_{\ell+1}^{(\mathrm{old})})^{\mathsf{T}}(BZ_{\ell_{0}}^{(\mathrm{new})}-Z_{\ell_{0}}^{(\mathrm{new})}\Theta^{(\mathrm{new})})=(Q_{\ell_{0}+1}^{(\mathrm{new})})^{\mathsf{T}}BZ_{\ell_{0}}^{(\mathrm{new})}$

is used.

###  4.2 Partial reorthogonalization

In floating-point arithmetic, round-off error compromises the orthogonality of
the basis computed by the block Lanczos method, leading to the emergence of
multiple converged Ritz values, so called ghost eigenvalues. To obtain
reliable information about the null space and its dimension, it is imperative
to maintain the orthogonality of Lanczos vectors with reorthogonalization
strategies.

Commonly employed reorthogonalization schemes include selective
reorthogonalization [24], partial reorthogonalization [31], and full
reorthogonalization. In our implementation, we use partial
reorthogonalization, which is cheaper than full reorthogonalization and more
effective than selective reorthogonalization [30]. Following the strategies
described in [33, 12], we define a $W$-recurrence to monitor loss of
orthogonality. For brevity, we omit the superscripts $\cdot^{(\mathrm{old})}$
and $\cdot^{(\mathrm{new})}$ regarding restarting. Additionally, we focus on
the reorthogonalization of the new basis vector $Q_{\ell+1}$ versus
$Z_{\ell}=[Q_{1},\dotsc,Q_{\ell}]$.

Supposing that $\ell_{0}d$ Ritz vectors were retained in the last restart, and

$\Theta=\mathrm{diag}\\{\Theta_{1},\dotsc,\Theta_{\ell_{0}}\\}\quad\text{and}\quad
F=[F_{1}^{\mathsf{T}},\dotsc,F_{\ell_{0}}^{\mathsf{T}}]^{\mathsf{T}},$

the recurrence Eq. 23 in finite precision can be written as

(25) $\displaystyle BQ_{i}$
$\displaystyle=Q_{i}\Theta_{i}+Q_{\ell_{0}}F_{i}^{\mathsf{T}}+\Delta_{i}\quad\text{for}\quad
1\leq i\leq\ell_{0},$ $\displaystyle BQ_{\ell_{0}+1}$
$\displaystyle=\sum_{i=1}^{\ell_{0}}Q_{i}F_{i}+Q_{\ell_{0}+1}\Xi_{\ell_{0}+1}+Q_{\ell_{0}+2}\Gamma_{\ell_{0}+1}^{\mathsf{T}}+\Delta_{\ell_{0}+1},$
$\displaystyle BQ_{j}$
$\displaystyle=Q_{j-1}\Gamma_{j-1}+Q_{j}\Xi_{j}+Q_{j+1}\Gamma_{j}^{\mathsf{T}}+\Delta_{j}\quad\text{for}\quad\ell_{0}+2\leq
j\leq\ell,$

where the matrices denoted with $\Delta$ account for round-off errors.
Following [12, 23], we may assume that the matrices $Q_{k}$ have orthonormal
columns and

(26) $\lVert\Delta_{k}\rVert\leq d\sqrt{n}\mathbf{u}\lVert
B\rVert=:\mathbf{u}_{d}$

for $1\leq k\leq\ell+1$, where $\mathbf{u}$ denotes the machine precision.
Different from full orthogonalization, we allow $Q_{i}$ and $Q_{j}$ to be not
(numerically) orthogonal to each other for $i\not=j$. The basic idea of
partial reorthogonalization is to monitor this loss of orthogonality by
examining $\lVert Q_{k}^{\mathsf{T}}Q_{\ell+1}\rVert$ for $1\leq k\leq\ell$.
If this value exceeds a certain threshold $\mathbf{u}_{\mathsf{RO}}$,
typically chosen as $\mathbf{u}_{\mathsf{RO}}=\sqrt{\mathbf{u}}$, we perform
reorthogonalization. Because computing $\lVert
Q_{k}^{\mathsf{T}}Q_{\ell+1}\rVert$ is computationally expensive, we instead
aim to calculate and use an upper bound for this value. Reorthogonalization is
not performed if the upper bound is below $\mathbf{u}_{\mathsf{RO}}$,
potentially cutting down computational costs while maintaining a decent level
of orthogonality.

To derive an upper bound for $\lVert Q_{k}^{\mathsf{T}}Q_{\ell+1}\rVert$, we
take $j=\ell$ in the third recurrence in Eq. 25. Multiplying with
$Q_{k}^{\mathsf{T}}$ from the left and rearranging terms yields the inequality

$\displaystyle\lVert Q_{k}^{\mathsf{T}}Q_{\ell+1}\rVert$
$\displaystyle\leq\lVert\Gamma_{\ell}^{-1}\rVert\bigl{\lVert}Q_{k}^{\mathsf{T}}BQ_{\ell}-Q_{k}^{\mathsf{T}}Q_{\ell-1}\Gamma_{\ell-1}-Q_{k}^{\mathsf{T}}Q_{\ell}\Xi_{\ell}-Q_{k}^{\mathsf{T}}\Delta_{\ell}\bigr{\rVert}$
$\displaystyle\leq\lVert\Gamma_{\ell}^{-1}\rVert\bigl{(}\lVert
Q_{k}^{\mathsf{T}}BQ_{\ell}\rVert+\lVert
Q_{k}^{\mathsf{T}}Q_{\ell-1}\rVert\lVert\Gamma_{\ell-1}\rVert+\lVert
Q_{k}^{\mathsf{T}}Q_{\ell}\rVert\lVert\Xi_{\ell}\rVert+\lVert
Q_{k}\rVert\lVert\Delta_{\ell}\rVert\bigr{)}.$

For $1\leq k\leq\ell_{0}$, the term $\lVert Q_{k}^{\mathsf{T}}BQ_{\ell}\rVert$
is bounded by using the symmetry of $B$ and the first recurrence in Eq. 25:

$\displaystyle\lVert Q_{k}^{\mathsf{T}}BQ_{\ell}\rVert$ $\displaystyle=\lVert
Q_{\ell}^{\mathsf{T}}BQ_{k}\rVert\leq\lVert
Q_{k}^{\mathsf{T}}Q_{\ell}\rVert\lVert\Theta_{k}\rVert+\lVert
Q_{\ell_{0}}^{\mathsf{T}}Q_{\ell}\rVert\lVert F_{k}\rVert+\lVert
Q_{\ell}\rVert\lVert\Delta_{k}\rVert.$

Combining these two inequalities and using Eq. 26 yields

$\lVert
Q_{k}^{\mathsf{T}}Q_{\ell+1}\rVert\leq\lVert\Gamma_{\ell}^{-1}\rVert\bigl{(}\lVert
Q_{k}^{\mathsf{T}}Q_{\ell}\rVert(\lVert\Theta_{k}\rVert+\lVert\Xi_{\ell}\rVert)+\lVert
Q_{\ell_{0}}^{\mathsf{T}}Q_{\ell}\rVert\lVert F_{k}\rVert+\lVert
Q_{k}^{\mathsf{T}}Q_{\ell-1}\rVert\lVert\Gamma_{\ell-1}\rVert+2\mathbf{u}_{d}\bigr{)}.$

We keep track of the recursive inequalities above using a strictly upper
triangular matrix $W$, such that each entry $w_{i,j}$ is an upper bound on
$\lVert Q_{i}^{\mathsf{T}}Q_{j}\rVert$ for $1\leq i<j\leq\ell+1$. Assuming
that columns $1,\ldots,\ell$ have been determined, the entries $1,\ldots,\ell$
of column $\ell+1$ are determined by the recurrences above as follows:

$w_{k,\ell+1}=\left\\{\begin{aligned}
\lVert\Gamma_{\ell}^{-1}\rVert\Bigl{(}&w_{\ell_{0}+1,\ell}(\lVert\Theta_{k}\rVert+\lVert\Xi_{\ell}\rVert)+w_{\ell_{0},\ell}\lVert
F_{k}\rVert\\\
&+w_{k,\ell-1}\lVert\Gamma_{\ell-1}\rVert+2\mathbf{u}_{d}\Bigr{)},&\\!\\!\\!\\!\\!\\!\\!\\!\text{for
}1\leq k\leq\ell_{0}.\\\
\lVert\Gamma_{\ell}^{-1}\rVert\Bigl{(}&w_{k,\ell}(\lVert\Xi_{k}\rVert+\lVert\Xi_{\ell}\rVert)+\sum_{1\leq
i\leq\ell_{0}}w_{i,\ell}\lVert F_{i}\rVert\\\
&+w_{k+1,\ell}\lVert\Gamma_{k}\rVert+w_{k,\ell-1}\lVert\Gamma_{\ell-1}\rVert+2\mathbf{u}_{d}\Bigr{)},&\\!\\!\\!\\!\\!\\!\text{for
}k=\ell_{0}+1.\\\
\lVert\Gamma_{\ell}^{-1}\rVert\Bigl{(}&w_{k,\ell}(\lVert\Xi_{k}\rVert+\lVert\Xi_{\ell}\rVert)+w_{k+1,\ell}\lVert\Gamma_{k}\rVert\\\
&+w_{k-1,\ell}\lVert\Gamma_{k-1}\rVert+w_{k,\ell-1}\lVert\Gamma_{\ell-1}\rVert+2\mathbf{u}_{d}\Bigr{)},&\text{for
}\ell_{0}+2\leq k\leq\ell.\end{aligned}\right.$

Note that the computation of the column $\ell+1$ only involves the previous
two columns $\ell-1$ and $\ell$, consistent with the three-term recurrence of
the Lanczos method. The entries of the initial two columns $\ell_{0}+1$ and
$\ell_{0}+2$ of $W$ are set to $\mathbf{u}_{d}$ because $Q_{\ell_{0}+1}$ and
$Q_{\ell_{0}+2}$ are reorthogonalized versus $Z_{\ell_{0}}$ immediately after
restarting. We monitor the loss of orthogonality using $W$. If
$w_{k,\ell+1}>\mathbf{u}_{\mathsf{RO}}$ for some $k$, we reorthogonalize
$Q_{\ell}$ and $Q_{\ell+1}$ versus $Z_{\ell-1}$ and set the strictly upper
triangular entries of the corresponding columns $\ell$ and $\ell+1$ of $W$ to
$\mathbf{u}_{d}$.

###  4.3 Pseudocode

Algorithm 2 incorporates the improvements to Algorithm 1 discussed in this
section. For simplicity, we assume that $\dim\mathcal{K}$ is an integer
multiple of $d$. The pseudocode uses Matlab’s colon notation to denote
submatrices.

Input: Matrix $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ and an unknown
number of zero singular values. Perturbation parameter $\epsilon>0$. Block
size $d$. Maximum allowed Krylov subspace dimension $\dim\mathcal{K}$.

Output: Orthonormal basis $V$ approximating the null space of $A$.

1 Set $B=A^{\mathsf{T}}A+\epsilon D$, with a diagonal matrix $D$ having
diagonal entries uniformly i.i.d. in $[0,1]$;

2 Draw Gaussian random matrix $\Omega\in\mathbb{R}^{n\times d}$;

3 Set $\mathbf{u}_{d}=d\sqrt{n}\mathbf{u}\lVert B\rVert$ and $\ell_{0}=0$;

4 Compute $Q_{1}$ as an orthonormal basis of $\Omega$;

5 while  _not converged_ do

6 Compute $Q_{\ell_{0}+1}^{\prime}=BQ_{\ell_{0}+1}$ and
$\Xi_{\ell_{0}+1}=Q_{\ell_{0}+1}^{\mathsf{T}}Q_{\ell_{0}+1}^{\prime}$;

7 Compute QR decomposition
$(I-Z_{\ell_{0}+1}Z_{\ell_{0}+1}^{\mathsf{T}})Q_{\ell_{0}+1}^{\prime}=Q_{\ell_{0}+2}\Gamma_{\ell_{0}+1}$;

8 Set $Z_{\ell_{0}+2}=[Z_{\ell_{0}+1},Q_{\ell_{0}+2}]$;

9 for  _$\ell=\ell_{0}+2,\dotsc,\dim\mathcal{K}/d$_ do

10 Compute $Q_{\ell}^{\prime}=BQ_{\ell}-Q_{\ell-1}\Gamma_{\ell-1}$ and
$\Xi_{\ell}=Q_{\ell}^{\mathsf{T}}Q_{\ell}^{\prime}$;

11 Compute QR decomposition
$Q_{\ell}^{\prime}-Q_{\ell}\Xi_{\ell}=Q_{\ell+1}\Gamma_{\ell}$;

12 Perform partial reorthogonalization as discussed in Section 4.2 and set
$Z_{\ell}=[Z_{\ell-1},Q_{\ell}]$ ;

13

14 end for

15 Form block tridiagonal matrix $T_{\ell}$ defined in Eq. 24 and compute its
eigenvalue decomposition $T_{\ell}V_{Z}=V_{Z}\Theta$ with eigenvalues in
ascending order;

16 Let $N_{\ell}$ be the number of eigenvalues of $T_{\ell}$ below $3\epsilon$
and set $\ell_{0}=(N_{\ell}+\dim\mathcal{K})/2$;

17 Set $\Theta=\Theta(1:\ell_{0},1:\ell_{0})$ and compute
$Z_{\ell_{0}}=Z_{\ell}V_{Z}(:,1:\ell_{0})$;

18 Compute QR decomposition
$(I-Z_{\ell_{0}}Z_{\ell_{0}}^{\mathsf{T}})Q_{\ell+1}=Q_{\ell_{0}+1}\Gamma_{\ell_{0}}$;

19 Set $F=([O,\Gamma_{\ell_{0}}]V_{Z})^{\mathsf{T}}$ and
$Z_{\ell_{0}+1}=[Z_{\ell_{0}},Q_{\ell_{0}+1}]$;

20

21 end while

22return _$V=Z_{\ell}V_{Z}(:,1:N_{\ell})$_ ;

Algorithm 2 Improved randomized small-block Lanczos

###  4.4 Preconditioning

In Remark 3.12, we established that the number of iterations $\ell$ needed for
convergence depends on the effective condition number
$\kappa=\overline{\sigma}/\underline{\sigma}$. Preconditioning can be used to
reduce this number.

####  4.4.1 Inner preconditioning

Given an $m\times m$ symmetric positive preconditioner $P_{L}$ for
$AA^{\mathsf{T}}$, one may replace $B=A^{\mathsf{T}}A+\epsilon D$ by

$B_{L}=A^{\mathsf{T}}P_{L}^{-1}A+\epsilon D$

within the Lanczos method. Besides maintaining the null space (the null spaces
of $A^{\mathsf{T}}P_{L}^{-1}A$ and $A$ are identical), the effective condition
number is improved to
$\kappa=\lVert(A^{\mathsf{T}}P_{L}^{-1}A)(A^{\mathsf{T}}P_{L}^{-1}A)^{\dagger}\rVert^{1/2}$.
For $m\gg n$, using such a preconditioner may be expensive because $P_{L}$ is
an $m\times m$ matrix.

####  4.4.2 Outer preconditioning

Given an invertible matrix $P_{R}\in\mathbb{R}^{n\times n}$ such that
$P_{R}P_{R}^{\mathsf{T}}$ is a preconditioner for $A^{\mathsf{T}}A$, one may
replace $B=A^{\mathsf{T}}A+\epsilon D$ by

$B_{R}=P_{R}^{-1}A^{\mathsf{T}}AP_{R}^{-\mathsf{T}}+\epsilon D$

within the Lanczos method. This transformation improves the effective
condition number to
$\kappa=\lVert(AP_{R}^{-\mathsf{T}})(AP_{R}^{-\mathsf{T}})^{\dagger}\rVert$.
However, unlike inner preconditioning, this method necessitates the
application of both $P_{R}^{-1}$ and $P_{R}^{-\mathsf{T}}$ in a single Lanczos
step. Moreover, an additional step of $P_{R}^{-\mathsf{T}}V$ is necessary to
post-process the output because $V$ approximates the null space of
$AP_{R}^{-\mathsf{T}}$ instead of $A$.

##  5 Numerical experiments

In this section, we present some numerical results to support the theoretical
results and demonstrate the efficiency of our method. All numerical
experiments111Scripts to reproduce numerical results are publicly available at
https://github.com/nShao678/Nullspace-code in this section are implemented in
Matlab 2022b and executed with an AMD Ryzen 9 6900HX Processor (8 cores,
3.3–4.9 GHz) and 32 GB of RAM.

###  5.1 Benchmark examples

For our experiments, we will consider two applications of null space
computation: (1) Counting the connected components of an undirected graph and
(2) cohomology computation.

For an undirected graph $G$, the null space $V_{G}$ of its graph Laplacian
$L_{G}$ contains information about its connectivity. Specifically, the
dimension of $V_{G}$ is the number of connected components, and the nodes in
each connected component can be computed from a null space basis. As a
specific example for a graph, we consider the mathematicians’ collaboration
network from MathSciNet [26]. As this graph is connected, the null space of
$L_{G}$ is trivial. To create a more interesting scenario, we consider the
sub-graph $G_{0}$ obtained by first deleting nodes of degree less than $N_{0}$
for a fixed positive integer $N_{0}$ (along with the corresponding edges) and
then deleting isolated nodes. Table 2 provides statistics in the resulting
graphs for $N_{0}=3$ and $N_{0}=10$. As the graph Laplacian $A$ is symmetric
positive semi-definite, we will replace $B=A^{\mathsf{T}}A+\epsilon D$ by
$B=A+\epsilon D$ in line 1 of Algorithm 2 for these examples. Moreover, in
order to accelerate convergence, we use outer preconditioning (see Section
4.4.2) with the preconditioner $P_{R}$ returned by the incomplete Cholesky
decomposition of $A$, using the Matlab function

(27)
$P_{R}=\mathtt{ichol}(A,\mathtt{struct}\text{(`type',`ict',`droptol',$10^{-3}$,`diagcomp',$0.1$)}).$

For $N_{0}=10$, this improves the effective condition number
$\kappa=\overline{\sigma}/\underline{\sigma}$ from $6227.1$ to $23.2$.

Table 2: Statistics on graph Laplacians for MathSciNet collaboration network after preprocessing. $N_{0}$ | #nodes | #connected components
---|---|---
3 | 168692 | 122
10 | 39890 | 66

For cohomology computations, we select the matrices from the SuiteSparse
collection [5] listed in Table 3.

Table 3: SuiteSparse collection matrices related to cohomology computations. Name | #rows | #columns | nnz | type | nullity
---|---|---|---|---|---
Franz7 | 10164 | 1740 | 40424 | Rings | 113
Franz8 | 16728 | 7176 | 100368 | Rings | 1713
Franz9 | 19588 | 4164 | 97508 | Rings | 619
Franz10 | 19588 | 4164 | 97508 | Rings | 619
GL7d12 | 8899 | 1019 | 37519 | $\mathrm{GL}_{7}(\mathbb{Z})$ | 59
GL7d13 | 47271 | 8899 | 356232 | $\mathrm{GL}_{7}(\mathbb{Z})$ | 961
GL7d14 | 171375 | 47271 | 1831183 | $\mathrm{GL}_{7}(\mathbb{Z})$ | 7939

###  5.2 Single-vector Lanczos

We first assess the convergence of the single-vector Lanczos method with full
reorthogonalization and without restarting for the graph Laplacian with
$N_{0}=10$ from Section 5.1, utilizing the incomplete Cholesky preconditioner
Eq. 27.

We varied the perturbation level $\epsilon$ from $4^{-4}$ to $4^{-10}$ and ran
single-vector Lanczos method for $\ell=1000$ iterations to ensure convergence.
Table 4 contains $\lambda_{N}(T_{\ell})$, the $N$th smallest eigenvalue of
$T_{\ell}$, and norms $\lVert V^{\mathsf{T}}AV\rVert$ and $\lVert AV\rVert$.
As expected, $\lambda_{N}(T_{\ell})$ and $\lVert V^{\mathsf{T}}AV\rVert$ are
observed to be of the order $\mathcal{O}(\epsilon)$ and
$\mathcal{O}(\epsilon^{2})$, respectively.

Table 4: Convergence of null space, where $\epsilon$ is taken as $4^{-k}$ for $k=4,\dotsc,10$. $\epsilon$ | 3.9e-03 | 9.8e-04 | 2.4e-04 | 6.1e-05 | 1.5e-05 | 3.8e-06 | 9.5e-07
---|---|---|---|---|---|---|---
$\lambda_{N}(T_{\ell})$ | 3.7e-03 | 9.2e-04 | 2.3e-04 | 5.8e-05 | 1.4e-05 | 3.6e-06 | 9.0e-07
$\lVert AV\rVert$ | 4.7e-03 | 1.2e-03 | 2.9e-04 | 7.4e-05 | 1.8e-05 | 4.6e-06 | 1.1e-06
$\lVert V^{\mathsf{T}}AV\rVert$ | 1.5e-06 | 9.2e-08 | 5.7e-09 | 3.6e-10 | 2.2e-11 | 1.4e-12 | 8.7e-14

The result of Theorem 3.10 predicts that the Ritz values of the single-vector
Lanczos method converge individually to the perturbed zero eigenvalues, that
is, some converge faster than others.. The numerical results shown in Fig. 1
for $\epsilon=4^{-10}$ confirm this result, showing that the number of
converged Ritz values grows steadily, a Ritz value converges approximately
every $10$ steps of the Lanczos method.

Figure 1: Convergence history of single-vector Lanczos with $\epsilon=4^{-10}$
for graph Laplacian: Number of converged Ritz values vs. number of iterations.

###  5.3 Small-block Lanczos method with restarting and partial
reorthogonalization

We now investigate the performance of Algorithm 2, the small-block Lanczos
method with restarting and partial reorthogonalization. For all remaining
experiments, we set the perturbation parameter to $\epsilon=10^{-3}$.

####  5.3.1 Comparison with null

First, we compare Algorithm 2 with the Matlab built-in function null, which is
based on the SVD, for the first six matrices from Table 3. The obtained
numerical results are presented in Table 5. For all parameter choices, the
computed nullity matches the actual nullity from Table 3. It can be observed
that restarting is very effective at limiting the dimension of the Krylov
subspace. Although increasing the block size increases the number of matrix-
vector multiplications, this increase is, for smaller $d$, outpaced by the
higher efficiency of block matrix-vector multiplications. Once the block size
becomes too large, both the iteration numbers and execution times start
increasing. Based on our experiments, a suitable strategy for choosing the
block size is $d=\dim\mathcal{K}/16$, where $\dim\mathcal{K}$ is the maximum
dimension of the Krylov subspace before restarting. Algorithm 2 is faster than
null for all examples, and the difference becomes more significant as the
matrix size increases.

Table 5: Number of matrix vector multiplications and execution time of
Algorithm 2 applied to cohomology computations. $\dim\mathcal{K}$ is the
maximum dimension of Krylov subspace before restarting. The last column refers
to Matlab’s built-in function null.

Block size | $d=2$ | $d=4$ | $d=8$ | $d=16$ | $d=32$ | null
---|---|---|---|---|---|---
$\\#$matvec | 2520 | 2616 | 2824 | 3232 | 3904 | $\times$
Time (s) | 0.9023 | 0.7106 | 0.5412 | 0.6125 | 0.6515 | 1.5449
(a)

Block size | $d=16$ | $d=32$ | $d=64$ | $d=128$ | $d=256$ | null
---|---|---|---|---|---|---
$\\#$matvec | 22848 | 22880 | 23808 | 25216 | 26880 | $\times$
Time (s) | 96.83 | 77.53 | 71.26 | 69.55 | 77.29 | 82.47
(b)

Block size | $d=8$ | $d=16$ | $d=32$ | $d=64$ | $d=128$ | null
---|---|---|---|---|---|---
$\\#$matvec | 12496 | 13008 | 12960 | 14528 | 16000 | $\times$
Time (s) | 17.68 | 13.63 | 11.11 | 9.852 | 9.572 | 21.47
(c)

Block size | $d=8$ | $d=16$ | $d=32$ | $d=64$ | $d=128$ | null
---|---|---|---|---|---|---
$\\#$matvec | 12440 | 12784 | 12960 | 14464 | 16000 | $\times$
Time (s) | 17.54 | 12.81 | 10.43 | 9.74 | 9.028 | 21.15
(d)

Block size | $d=1$ | $d=2$ | $d=4$ | $d=8$ | $d=16$ | null
---|---|---|---|---|---|---
$\\#$matvec | 1439 | 1482 | 1504 | 1720 | 1920 | $\times$
Time (s) | 0.2188 | 0.1858 | 0.1545 | 0.1752 | 0.2387 | 0.6476
(e)

Block size | $d=16$ | $d=32$ | $d=64$ | $d=128$ | $d=256$ | null
---|---|---|---|---|---|---
$\\#$matvec | 17552 | 18176 | 19072 | 19712 | 23040 | $\times$
Time (s) | 55.20 | 44.71 | 38.81 | 37.29 | 41.50 | 212.2
(f)

####  5.3.2 Comparison with large-block Lanczos

Traditionally, the block Lanczos method with restarting is used with a block
size that is slightly larger than the nullity. Obviously, when memory is a
constraint, the small-block or even single-vector Lanczos methods with
restarting are preferred. Even when memory is not a constraint, using the
small-block version can be beneficial. To illustrate this, we consider two
examples: the matrix “Franz9” from cohomology computations and the graph
Laplacian with $N_{0}=3$, which have nullities $619$ and $122$, respectively.
We compare the performance of the small-block and large-block Lanczos methods.
For the large-block Lanczos method, we set the maximum dimension of the Krylov
subspace before restarting to be approximately four times the nullity, i.e.,
$2560$ and $512$, respectively.

As shown in Table 6, the large-block Lanczos method requires more matrix-
vector multiplications and execution time compared to the small-block version.
The nullity has been computed correctly by all methods. One advantage of the
large-block Lanczos method is that it does not require small perturbations
such as $B=A^{\mathsf{T}}A+\epsilon D$ or $B=A+\epsilon D$, thanks to the use
of a sufficiently large-block size. While the absence of perturbations yields
more accurate results, the overall computational cost remains unaffected
because the additional computational cost from applying the diagonal
perturbation matrix $D$ is negligible.

Table 6: Number of matrix vector multiplications, execution time, and null
space residual of Algorithm 2 applied to “Franz9” and the graph Laplacian with
$N_{0}=3$. $\dim\mathcal{K}$ is the maximum dimension of Krylov subspace
before restarting. The last column refers to the block Lanczos method applied
to the original matrix (without diagonal perturbations). The algorithm is
stopped once $\lVert AV\rVert\leq 10^{-4}$ for “Franz9” or $\lVert
AV\rVert\leq 1.2\times 10^{-3}$ for the graph Laplacian is detected before
restarting.

| Small-block Lanczos | Block Lanczos
---|---|---
Block size | $d=40$ | $d=80$ | $d=160$ | $d=320$ | $d=640$ | $d=640$
$\\#$matvec | 7960 | 6960 | 7200 | 8960 | 10240 | 10240
Time (s) | 13.57 | 10.27 | 9.925 | 12.32 | 11.23 | 11.38
$\lVert AV\rVert$ | 9.6e-05 | 9.6e-05 | 1.0e-04 | 9.6e-05 | 9.6e-05 | 1.2e-05
(g)

| Small-block Lanczos | Block Lanczos
---|---|---
Block size | $d=8$ | $d=16$ | $d=32$ | $d=64$ | $d=128$ | $d=128$
$\\#$matvec | 1608 | 1632 | 1664 | 2048 | 3584 | 3584
Time (s) | 48.55 | 43.52 | 40.29 | 44.42 | 70.72 | 69.70
$\lVert AV\rVert$ | 1.0e-03 | 1.0e-03 | 1.0e-03 | 1.1e-03 | 1.0e-03 | 7.5e-04
(h)

####  5.3.3 Large-scale problems

In the final experiment, we address a particularly challenging problem where
both the dense method and the traditional block Lanczos method are infeasible
due to memory constraints. For the matrix “GL7d14”, it is impossible to apply
the null function because the dense matrix cannot be stored in memory.
Similarly, the large-block Lanczos method encounters a memory bottleneck,
making it impractical to perform a restart after even just two iterations.
Table 7 demonstrates that the the small-block Lanczos method successfully
computes the null space within a Krylov subspace of dimension 10240.
Remarkably, this dimension is only 1.29 times larger than the nullity,
highlighting the efficiency and feasibility of the small-block method in
handling such large-scale problems.

Table 7: Number of matrix vector multiplications and execution time of Algorithm 2 applied to “GL7d14” with nullity $7939$. The maximum dimension of Krylov subspace before restarting is $10240$. The computed nullity is correct and $\lVert AV\rVert\approx 1.4\times 10^{-4}$ Block size | $d=16$ | $d=32$ | $d=64$ | $d=128$ | $d=256$ | $d=512$
---|---|---|---|---|---|---
$\\#$matvec | 133456 | 133344 | 132288 | 132352 | 133888 | 140800
Time(s) | 11516 | 8137.4 | 6986.4 | 6375.0 | 6218.7 | 6716.4

##  6 Conclusion

In this work, we proposed a novel randomized small-block Lanczos method that
incorporates modern Krylov subspace techniques such as restarting and partial
reorthogonalization to compute null spaces of large sparse matrices. For the
examples considered, our implementation outperforms both the default SVD-based
null space solver and the traditional block Lanczos method. Additionally, we
provide a theoretical analysis of the convergence for the single-vector case,
that is, $d=1$. It remains an open problem to derive convergence bounds that
reflect the good performance observed for small $d>1$.

## Acknowledgments

The authors thank Yuji Nakatsukasa and Rikhav Shah for helpful discussions on
this work.

## References

  * [1] M. Aizenman, R. Peled, J. Schenker, M. Shamis, and S. Sodin, Matrix regularizing effects of Gaussian perturbations, Communications in Contemporary Mathematics, 19 (2017), p. 1750028, https://doi.org/10.1142/S0219199717500286.
  * [2] K. Batselier, P. Dreesen, and B. De Moor, On the null spaces of the Macaulay matrix, Linear Algebra and its Applications, 460 (2014), pp. 259–289, https://doi.org/https://doi.org/10.1016/j.laa.2014.07.035.
  * [3] M. Berry, M. Heath, I. Kaneko, M. Lawo, R. Plemmons, and R. Ward, An algorithm to compute a sparse basis of the null space, Numerische Mathematik, 47 (1985), pp. 483–504, https://doi.org/10.1007/BF01389453.
  * [4] T. A. Davis, Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization, ACM Transactions on Mathematical Software, 38 (2011), https://doi.org/10.1145/2049662.2049670.
  * [5] T. A. Davis and Y. Hu, The University of Florida sparse matrix collection, ACM Transactions on Mathematical Software (TOMS), 38 (2011), pp. 1–25, https://doi.org/10.1145/2049662.2049663.
  * [6] P. Drineas and I. C. Ipsen, Low-rank matrix approximations do not need a singular value gap, SIAM Journal on Matrix Analysis and Applications, 40 (2019), pp. 299–319, https://doi.org/10.1137/18M1163658.
  * [7] H. Edelsbrunner and J. L. Harer, Computational topology: an introduction, American Mathematical Society, Rhode Island, 2022.
  * [8] L. V. Foster, Rank and null space calculations using matrix decomposition without column interchanges, Linear Algebra and its Applications, 74 (1986), pp. 47–71, https://doi.org/10.1016/0024-3795(86)90115-1.
  * [9] L. V. Foster and T. A. Davis, Algorithm 933: Reliable calculation of numerical rank, null space bases, pseudoinverse solutions, and basic solutions using suitesparseQR, ACM Transactions on Mathematical Software, 40 (2013), https://doi.org/10.1145/2513109.2513116.
  * [10] G. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Maryland, 4th ed., 2013.
  * [11] C. Gotsman and S. Toledo, On the computation of null spaces of sparse rectangular matrices, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 445–463, https://doi.org/10.1137/050638369.
  * [12] R. G. Grimes, J. G. Lewis, and H. Simon, A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems, SIAM Journal on Matrix Analysis and Applications, 15 (1994), pp. 228–272, https://doi.org/10.1137/S0895479888151111.
  * [13] R. J. Hewett, M. T. Heath, M. D. Butala, and F. Kamalabadi, A robust null space method for linear equality constrained state estimation, IEEE Transactions on Signal Processing, 58 (2010), pp. 3961–3971, https://doi.org/10.1109/TSP.2010.2048901.
  * [14] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, New York, 2 ed., 2012, https://doi.org/10.1017/9781139020411.
  * [15] B. Laurent and P. Massart, Adaptive estimation of a quadratic functional by model selection, The Annals of Statistics, 28 (2000), pp. 1302 – 1338, https://doi.org/10.1214/aos/1015957395.
  * [16] R.-C. Li and L.-H. Zhang, Convergence of the block Lanczos method for eigenvalue clusters, Numerische Mathematik, 131 (2015), pp. 83–113, https://doi.org/10.1007/s00211-014-0681-6.
  * [17] R. Mathias, Quadratic residual bounds for the Hermitian eigenvalue problem, SIAM J. Matrix Anal. Appl., 19 (1998), pp. 541–550, https://doi.org/10.1137/S0895479896310536.
  * [18] R. Merris, Laplacian matrices of graphs: a survey, Linear Algebra and its Applications, 197 (1994), pp. 143–176, https://doi.org/10.1016/0024-3795(94)90486-3.
  * [19] R. Meyer, C. Musco, and C. Musco, On the unreasonable effectiveness of single vector krylov methods for low-rank approximation, in Proceedings of the 2024 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2024, pp. 811–845, https://doi.org/10.1137/1.9781611977912.32.
  * [20] N. Minami, Local fluctuation of the spectrum of a multidimensional Anderson tight binding model, Communications in Mathematical Physics, 177 (1996), pp. 709–725, https://doi.org/10.1007/BF02099544.
  * [21] C. Musco and C. Musco, Randomized block Krylov methods for stronger and faster approximate singular value decomposition, Advances in neural information processing systems, 28 (2015), https://proceedings.neurips.cc/paper_files/paper/2015/file/1efa39bcaec6f3900149160693694536-Paper.pdf.
  * [22] T. Park and Y. Nakatsukasa, A fast randomized algorithm for computing an approximate null space, BIT Numerical Mathematics, 63 (2023), p. 36, https://doi.org/10.1007/s10543-023-00979-7.
  * [23] B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, Philadelphia, 1998, https://doi.org/10.1137/1.9781611971163.
  * [24] B. N. Parlett and D. S. Scott, The Lanczos algorithm with selective orthogonalization, Mathematics of Computation, 33 (1979), pp. 217–238, https://doi.org/10.1090/S0025-5718-1979-0514820-3.
  * [25] R. Platt, A. H. Fagg, and R. A. Grupen, Null-space grasp control: Theory and experiments, IEEE Transactions on Robotics, 26 (2010), pp. 282–295, https://doi.org/10.1109/TRO.2010.2042754.
  * [26] R. A. Rossi and N. K. Ahmed, The network data repository with interactive graph analytics and visualization, in AAAI, 2015, https://networkrepository.com.
  * [27] Y. Saad, On the rates of convergence of the Lanczos and the block-Lanczos methods, SIAM Journal on Numerical Analysis, 17 (1980), pp. 687–706, https://doi.org/10.1137/0717059.
  * [28] Y. Saad, Numerical methods for large eigenvalue problems: revised edition, SIAM, Philadelphia, 2011, https://doi.org/10.1137/1.9781611970739.
  * [29] G. Shklarski and S. Toledo, Computing the null space of finite element problems, Computer Methods in Applied Mechanics and Engineering, 198 (2009), pp. 3084–3095, https://doi.org/10.1016/j.cma.2009.05.012.
  * [30] H. Simon, The Lanczos algorithm for solving symmetric linear systems, PhD thesis, University of California, Berkeley, 1982.
  * [31] H. Simon, The Lanczos algorithm with partial reorthogonalization, Mathematics of Computation, 42 (1984), pp. 115–142, https://doi.org/10.1090/S0025-5718-1984-0725988-X.
  * [32] G. W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM Journal on Matrix Analysis and Applications, 23 (2002), pp. 601–614, https://doi.org/10.1137/S0895479800371529.
  * [33] K. Wu and H. Simon, Thick-restart Lanczos method for large symmetric eigenvalue problems, SIAM Journal on Matrix Analysis and Applications, 22 (2000), pp. 602–616, https://doi.org/10.1137/S0895479898334605.
  * [34] J. Ye, T. Xiong, and D. Madigan, Computational and theoretical analysis of null space and orthogonal linear discriminant analysis., Journal of Machine Learning Research, 7 (2006), http://jmlr.org/papers/v7/ye06a.html.
  * [35] Q. Yuan, M. Gu, and B. Li, Superlinear convergence of randomized block lanczos algorithm, in 2018 IEEE International Conference on Data Mining (ICDM), IEEE, 2018, pp. 1404–1409, https://doi.org/10.1109/ICDM.2018.00193.
  * [36] L. Zhang, T. Xiang, and S. Gong, Learning a discriminative null space for person re-identification, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016, https://openaccess.thecvf.com/content_cvpr_2016/html/Zhang_Learning_a_Discriminative_CVPR_2016_paper.html.
  * [37] Y. Zhou and Y. Saad, Block Krylov-Schur method for large symmetric eigenvalue problems, Numerical Algorithms, 47 (2008), pp. 341–359, https://doi.org/10.1007/s11075-008-9192-9.
  * [38] P. Zhu and A. Knyazev, Angles between subspaces and their tangents, Journal of Numerical Mathematics, 21 (2013), pp. 325–340, https://doi.org/10.1515/jnum-2013-0013.