# Rigorous screened interactions for realistic correlated electron systems

Charles J. C. Scott<EMAIL_ADDRESS>Department of Physics, King’s
College London, Strand, London WC2R 2LS, United Kingdom  George H. Booth
<EMAIL_ADDRESS>Department of Physics, King’s College London, Strand,
London WC2R 2LS, United Kingdom

###### Abstract

We derive a widely-applicable first principles approach for determining two-
body, static effective interactions for low-energy Hamiltonians with
quantitative accuracy. The algebraic construction rigorously conserves all
instantaneous two-point correlation functions in a chosen model space at the
level of the random phase approximation, improving upon the traditional
uncontrolled static approximations. Applied to screened interactions within a
quantum embedding framework, we demonstrate these faithfully describe the
relaxation of local subspaces via downfolding high-energy physics in molecular
systems, as well as enabling a systematically improvable description of the
long-range plasmonic contributions in extended graphene.

Building effective models that retain the physics of interest but strip away
extraneous complexity is central to progress in understanding physical
mechanisms and emergent behaviour in complex systems. Nowhere is this more
true than interacting electron systems, where effective Hamiltonians occupy a
central position in both condensed matter and quantum chemistry, from the
Hubbard model to ligand field theory [1, 2, 3, 4]. A push is now underway to
erode the boundary between phenomenological effective models with empirical
parameters, and ab initio modelling with material specificity, where the act
of constructing a material-specific effective Hamiltonian is increasingly the
first step in a larger workflow involving its subsequent solution within a
multi-method approach. This is critical to extend the scope of accurate yet
computationally demanding many-body methods, placing significant urgency on
approaches in which the relevant physics outside a low-energy model space is
rigorously downfolded or renormalized into the effective Hamiltonian [5, 6,
7]. In this work we propose a simple and efficient approach to effective
interactions, demonstrating quantitative accuracy and specificity resulting
from integrating out ab initio high-energy and long-range physics, motivated
by the exact conservation of instantaneous expectation values.

A salient example in the need for accurate renormalized interactions is
quantum embedding, describing correlated many-body phenomena within a local
subspace [8, 9, 10, 11]. The missing interactions with states outside the
subspace should renormalize the effective subspace interactions in real
materials, generally necessitating the use of effective interactions that are
often parameterized empirically via Hubbard or Hund terms [12, 13, 14], or
downfolded from other theories [11, 15, 16, 17]. In quantum chemistry,
analogous subspaces are often described by a ‘complete active space’ (CAS),
whereby a small number of low-energy mean-field orbitals are chosen for an
accurate treatment of the strong correlation in ‘multi-reference’ approaches
[18]. More broadly, a wide range of both qualitative and quantitative studies
into correlated materials rely on first obtaining appropriately screened
effective interactions of a simplified model from an ab initio starting point
[19, 5].

It is often argued that the random phase approximation (RPA) contains much of
the physics missing from low-energy models, and is the appropriate theory in
which to construct effective interactions [20, 11, 13, 21]. As an infinite
resummation of the bubble diagrams, it correctly describes the long-range
charge fluctuations, plasmons, high-energy scattering and many-body dispersion
that predominates in high-density metallic and semi-conducting extended
systems [22, 23, 24]. While its traditional formulation lacks exchange or
ladder diagrams, their contributions generally decay more rapidly and can
often be captured within the model space, with approaches to screen beyond RPA
still an active research area [25, 15, 16].

The constrained RPA (cRPA) method has therefore become a widespread choice for
deriving low-energy effective interactions from first principles, applied from
molecules to Mott insulators and routinely as a precursor to quantum embedding
methods [26, 27, 13, 25, 28, 16, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
8, 40, 5, 41]. Since the RPA is a well-defined diagrammatic theory, the
bubbles corresponding to the polarizability within the model space
($\Pi_{m}(\omega)$) can be removed from the total polarization ($\Pi_{\rm
ext}(\omega)=\Pi(\omega)-\Pi_{m}(\omega)$). The screened interactions,
$\mathcal{U}(\omega)$, are then found from the infinite RPA resummation
coupling this space to the external degrees of freedom via this scattering
channel, as

$\displaystyle\mathcal{U}(\omega)=v+v\Pi_{\rm ext}(\omega)v+v\Pi_{\rm
ext}(\omega)v\Pi_{\rm ext}(\omega)v+\dots,$ (1)

where $v$ denotes the bare Coulomb interaction. In this way, double counting
of correlated effects are avoided once the resulting effective model is
solved, yet direct scattering events to all orders between the two spaces are
included [26]. The cRPA approach is also internally consistent: an RPA
calculation on the resulting effective Hamiltonian with interactions
$\mathcal{U}(\omega)$, will give exactly the results of the RPA on the full
system. More specifically, the cRPA ensures that the RPA density-density (dd)
response function, $\chi(\omega)$, of the model space with effective
interactions $\mathcal{U}(\omega)$, is identical to the projection of the RPA
full system dd-response into the model space. In this way, the effective
interactions can also be seen to be ‘state universal’, providing the correctly
renormalized interactions for the entire RPA spectrum of the model subspace.

While this puts the cRPA approach on a solid footing, few methods are
computationally tractable with the resulting dynamical interactions of Eq. 1.
Practical necessity therefore forces the widespread approximation of taking
the static, $\omega\rightarrow 0$ limit of the screened interaction of Eq. 1
(static-cRPA). This uncontrolled approximation is qualitatively justified in
capturing the relevant long-wavelength behaviour where the RPA is accurate
[25], however forces us to entirely give up on rigorous conservation of any
expectation values from the RPA. There is significant scepticism of its
accuracy, with the cRPA interactions often over-screening the physics [25, 28,
42, 15, 16]. Far from purely a small quantitative shift, this can result in
wrong energy ordering of states and phases in the resulting model, or missing
spectral weight transfer to plasmon satellites [43, 14, 35].

A further serious technical limitation of the cRPA approach is that the model
space must be chosen from a selection of low-energy bands, rather than
directly as local degrees of freedom. Specifically, it requires the
irreducible polarizability (defined by the reference mean-field) to have no
coupling between excitations within the model space and ones in the rest of
the system, otherwise $\mathcal{U}(\omega)$ alone is not sufficient to
reproduce $\chi(\omega)$. Effective local interactions are therefore found via
localization after determining the screening of these low-energy bands.
However, this is problematic, missing out on screening within the low-energy
bands themselves if the resulting local interactions are subsequently
truncated (e.g. to Hubbard/Hund form), and causing further difficulties when
it is not possible to fully represent the relevant subspace (e.g. atomic
$d$-shells) in the low energy bands due to significant hybridization with
other states[44, 29]. In addition, the resulting cRPA screened interactions
modify the reference mean-field state, resulting in unintended changes to the
subspace density and bandstructure [28, 42, 45].

To motivate our approach, we ask an alternative pertinent question: What
physical quantities can we rigorously conserve at the RPA level in a chosen
model space, under the constraint that the resulting effective interactions
must remain static and two-body? This leads to the development of the ‘moment-
constrained’ RPA (mRPA) approach, which exactly conserves the instantaneous
part of the two-point dd-response in the model subspace, as well as the
reference state. We argue that conserving physical expectation values in the
construction of effective static interactions provides a more rigorous
foundation than the widespread static-cRPA approximations.

_Moment-constrained RPA:-_ The RPA can be formulated as a quasi-bosonic
eigenvalue problem in the space of particle-hole excitations and de-
excitations of a reference state [46],

$\displaystyle\begin{bmatrix}\mathbf{A}&\mathbf{B}&\\\
\mathbf{-B}&\mathbf{-A}&\end{bmatrix}\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}=\begin{bmatrix}\mathbf{X}&\mathbf{Y}\\\
\mathbf{Y}&\mathbf{X}\end{bmatrix}\begin{bmatrix}\mathbf{\Omega}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{\Omega}\end{bmatrix},$ (2)

with all blocks being of dimension given by the product of the number of hole
and particle states. We define $\mathbf{A}=\mathbf{\Delta}+v$, where
$\mathbf{\Delta}$ is a diagonal matrix of particle-hole excitation energies
(provided by a mean-field), defining the poles of $\Pi(\omega)$, with
$\mathbf{B}=v$ providing the coupling between the excitations and de-
excitations via the bare Coulomb interaction, $v$, in the particle-hole ($ph$)
channel. The diagonal matrix $\mathbf{\Omega}$ provides the poles of the RPA
dd-response function, $\chi(\omega)$, with the residues defined by the
amplitudes of the excitations and de-excitations, $\mathbf{X}$ and
$\mathbf{Y}$ respectively. These form a biorthogonal set of eigenvectors,
providing the relations
$(\mathbf{X}+\mathbf{Y})^{-1}=(\mathbf{X}-\mathbf{Y})^{T}$ and
$(\mathbf{X}-\mathbf{Y})^{-1}=(\mathbf{X}+\mathbf{Y})^{T}$. We can expand
$\chi(\omega)$ as a Laurent series, with its dynamics fully characterized by
the moments of its spectral distribution [40, 47], as

$\eta^{(n)}=(\mathbf{X}+\mathbf{Y})\mathbf{\Omega}^{n}(\mathbf{X}+\mathbf{Y})^{T}\quad;\quad
n\in{\mathbb{Z}}.$ (3)

The zeroth moment of the distribution, $\eta^{(0)}$, characterizes the
instantaneous part of the correlated dd-response, as
$\langle({\hat{c}}^{\dagger}_{i}{\hat{c}}_{a}+{\hat{c}}_{a}^{\dagger}{\hat{c}}_{i})({\hat{c}}_{j}^{\dagger}{\hat{c}}_{b}+{\hat{c}}_{b}^{\dagger}{\hat{c}}_{j})\rangle-\langle{\hat{c}}^{\dagger}_{i}{\hat{c}}_{a}+{\hat{c}}_{a}^{\dagger}{\hat{c}}_{i}\rangle\langle{\hat{c}}_{j}^{\dagger}{\hat{c}}_{b}+{\hat{c}}_{b}^{\dagger}{\hat{c}}_{j}\rangle$,
summed over the same-spin particle-hole (de)excitations, denoted by indices
$(a,i)$ and $(b,j)$. This describes the correlated contribution to the two-
body reduced density matrix, and all resulting static expectation values at
the RPA level [48, 49, 50]. It is this quantity (through to first order) which
we aim to rigorously conserve within a chosen model space with our effective
interactions. Crucially, this can be achieved while maintaining static
renormalized effective interactions in the model space, preserving symmetries,
and without changing the reference mean-field in the model space.

The structure of the RPA equations imposes a relation between the first two
dd-moments [47],

$\displaystyle\eta^{(1)}$
$\displaystyle=(\mathbf{A}-\mathbf{B})=\eta^{(0)}(\mathbf{X}-\mathbf{Y})\mathbf{\Omega}(\mathbf{X}-\mathbf{Y})^{T}\eta^{(0)}$
(4) $\displaystyle=\eta^{(0)}(\mathbf{A}+\mathbf{B})\eta^{(0)}.$ (5)

Inserting the definitions of $\mathbf{A}$ and $\mathbf{B}$, we find an
equation linear in the interaction. This can be analytically inverted to find
an interaction kernel that under the RPA gives rise to a desired $\eta^{(0)}$
and $\eta^{(1)}$,

$v=\frac{1}{2}\left((\eta^{(0)})^{-1}\eta^{(1)}(\eta^{(0)})^{-1}-\mathbf{\Delta}\right).$
(6)

Therefore, by substituting all quantities on the RHS of Eq. 6 for their
projection of the full RPA into the chosen model space, an effective static
model space interaction can be found, $\mathcal{U}_{\textrm{mRPA}}$. This
ensures that the model space RPA with
$v\rightarrow\mathcal{U}_{\textrm{mRPA}}$ rigorously conserves all full system
RPA expectation values derived from $\mathbf{\eta}^{(0)}$ and
$\mathbf{\eta}^{(1)}$ by construction, which includes all instantaneous
correlators in the model space. A similar approach to conserve higher-order
moments results in equations non-linear in the resulting interaction, and
subsequently necessitates a dynamical component to the resulting interaction,
as expected to recover the dynamical cRPA limit of the full dd-response [40].
Focusing on conservation of just the first two dd-moments enables a fully
static model space screened interaction.

The model space with $\mathcal{U}_{\textrm{mRPA}}$ can also reproduce the
subspace contribution to the RPA correlation energy due to its dependence on
$\mathbf{\eta}^{(0)}$, as
$E_{\mathrm{corr}}^{\mathrm{RPA}}=\frac{1}{2}\mathrm{Tr}[\mathbf{\eta}^{(0)}(\mathbf{A}+\mathbf{B})-\mathbf{A}]$,
exploited in App. D to compute energetic corrections to the model space.
Beyond this, the first two spectral moments of the $GW$ self-energy in the
subspace are also exactly described with mRPA interactions, indicating a
formal conservation of certain correlated one-body properties [47].
Symmetries, including spin-independence of the resulting effective
interactions are also exactly preserved, derived in App. B. These rigorous
properties provide a robust footing for use of $\mathcal{U}_{\textrm{mRPA}}$
in subsequent correlated treatments.

Building $\mathcal{U}_{\textrm{mRPA}}$ via Eq. 6 requires the component of
$\mathbf{\eta}^{(0)}$ and $\mathbf{\Delta}$ in the subspace. Construction
directly from Eqs. 2-3 entails a prohibitive $\mathcal{O}\left(N^{6}\right)$
scaling, but this is reduced to an efficient $\mathcal{O}\left(N^{4}\right)$
following the approach of App. A. The only constraint on the choice of
subspace that the cluster excitation space is an orthogonal projection of the
full system excitation space, i.e. that particle- and hole-like character of
the subspace orbitals as defined by the reference state is preserved. While
this is trivially true where cRPA is valid via a subspace selected from mean-
field bands, this is a far looser requirement allowing for a subspace where a
non-diagonal $\Pi(\omega)$ couples the subspace to its environment to be
considered. This admits direct mRPA screening of arbitrary (e.g. local atomic-
like) subspaces by at most doubling their size by conserving the reference
mean-field density matrix over the subspace[51, 52], or local spaces formed by
localizing hole and virtual bands separately before screening. This direct
screening of local subspaces therefore includes screening via long-range low-
energy band transitions precluded in traditional cRPA, enabling direct
application to ab initio quantum embedding clusters.

Interestingly, the resulting $\mathcal{U}_{\textrm{mRPA}}$ only screens
interactions in the $ph$ channel (which is expected to be the dominant long-
range contribution to screening), rather than the full four-point interaction.
Other approximations (e.g. $T$-matrix for $pp$ diagrams [53, 54, 55, 56])
would screen other interaction channels in an analogous formulation, and
future work can consider the effect on mRPA interactions from these other
channels. It is justified that only the $ph$ channel interaction is screened
in mRPA unlike cRPA, since RPA itself is fully determined by this component of
the interaction, and mRPA is constructed to describe this RPA physics as
opposed to the use of the screening equation (Eq. 1) in cRPA. This feature of
only screening the $ph$ interactions, along with the conservation of the
reference density, ensures that the occupied reference bandstructure
(Hartree–Fock or Kohn–Sham) is unchanged with mRPA screened interactions.

Figure 1: Trace of the screened and bare interaction, $\mathcal{U}_{pq,pq}$ in
a (10,10) CAS space of Benzene in a cc-pVDZ basis. Shown are the dynamical
cRPA interactions, as well as the widely used static limit, and the static-by-
construction mRPA interactions. Inset: Error in the RPA moments as an
expansion of the dynamical subspace dd-response. The first two mRPA dd-moments
are exact by construction, yet also exhibit marginally smaller errors for
higher moments compared to the static cRPA interaction. Moment error is
computed as the mean squared error of $\mathbf{\eta}^{(n)}$ over the subspace,
normalized by the exact $||\mathbf{\eta}^{(n)}||_{2}$ at each order.

_Results:-_ We begin by benchmarking the screening of the low-energy
Hartree–Fock orbitals of Benzene (an active space of five hole and five
particle states). Benzene has previously been considered a paradigmatic
example in this context [25], while also small enough to enable comparison to
high-level reference results. In Fig. 1 we show the bare Coulomb, dynamical
and static-cRPA, and mRPA subspace interactions traced over all channels. The
$\mathcal{U}_{\textrm{mRPA}}$ is generally less screened than static-cRPA,
noting the recent evidence that static-cRPA overscreens interactions [28, 42,
15]. We can consider the accuracy of resulting subspace observables by
comparing to coupled-cluster (CCSD)[57, 58, 59]. CCSD encodes the wavefunction
in $T$-amplitudes, and we can consider the projection of the full-system
$T$-amplitudes into the active space as an accurate subspace description. This
is used to benchmark subspace-only CCSD calculations with the different static
interactions 111Coupled-cluster is defined in terms of the Fock matrix and
fluctuation potential. The virtual-virtual block of the Fock matrix is shifted
by $\Delta f_{ab}=\sum_{i\in\textrm{occ}}(ia|bi)-\mathcal{U}_{ia,ib}$ due to
the screened $ph$ mRPA interactions, but other blocks are unchanged ensuring
the reference-independence of mRPA.. We find the mean squared error of the
subspace CCSD $T_{2}$-amplitudes to be 1.10 for bare subspace interactions,
reduced to 0.25 for static-cRPA and only 0.20 for mRPA, indicating that all
ground-state expectation values within the subspace can expected to be more
faithfully reproduced with mRPA interactions than the alternatives.

While we expect ground states to be particularly faithful given the mRPA
conservation of instantaneous correlators, we can also compare the ability of
these effective interactions to reproduce the full subspace excitation
spectrum. Fig. 1 quantifies this via the errors in the subspace RPA moment
expansion of Eq. 3 which fully characterizes the dd-response. While the first
two moments are exact for mRPA by construction (and $\mathbf{\eta}^{(1)}$ for
static-cRPA), the errors in higher moments are also marginally reduced in
their relative error compared to static-cRPA, indicating that the fidelity of
the subspace dd-response over all frequencies is at least as accurate as
static-cRPA.

We now consider the diverse W4-11 test set of 150 molecules, exhibiting a wide
range of bonding, radical and correlated physics [61]. Moving towards direct
application of mRPA to quantum embedding methodologies, we consider
fragmenting each molecule into individual atoms comprising a minimal set of
their intrinsic valence atomic orbitals (IAOs)[62]. These IAO fragments are
augmented with an interacting bath of at most the same dimension as the IAO
fragment, following the static density matrix embedding (DMET) approach[51,
63, 52, 64, 65, 66]. This defines multiple small atomic subspaces for each
molecule, which can be individually solved via CCSD with either bare or mRPA
interactions. The resulting subspace CCSD states are recombined to provide a
total energy estimator over the subspaces for each molecule (see Refs. 66, 67
for more details). These are compared to the energy of the exact full-system
CCSD projected into these subspaces. This discrepancy quantifies the ability
of the mRPA to account for relaxation of these embedded atomic descriptions
due to the neglected interactions with the rest of the molecule.

Figure 2: Distribution of energy errors in atomically fragmented DMET-CCSD
with bare or mRPA interactions across the W4-11 test set of 150 molecules, in
increasingly large basis sets. Errors are computed comparing to the full-
system CCSD of each molecule projected into each subspace.

Figure 2 shows this error aggregated over the test set in increasingly large
basis sets, where screening of these low-energy subspaces by higher-energy
scattering dominates. The mRPA interactions screen the fragments defined by
the embedding, resulting in a substantial $\sim 66\%$ reduction in both the
mean absolute error and standard deviation of the energy error across the data
set. This consistent reduction in error across larger basis sets points for
these diverse systems attests to the broad applicability of the mRPA
interactions for correlated systems. We note that cRPA cannot be easily
compared in this context, due to the difficulties in directly screening
subspaces that are not mean-field orbitals, as discussed previously.

Finally, we consider the application of mRPA to extended systems where long-
range collective plasmons strongly renormalize local properties, which are
difficult to describe in a local subspace[68, 25]. Furthermore, we demonstrate
systematic improvability of this screened local subspace. This is achieved in
a consistent framework via extending the bath space, formally including
additional states which exactly minimize the error of the subspace RPA
$\mathbf{\eta}^{(0)}$ with bare interactions, thereby spanning physics at
longer length-scales in the model subspaces. The algebraic construction of
these additional bath states relies on evaluation of the same quantities as
the mRPA interactions, and systematically enlarges the local correlated
subspaces in an optimal way to completeness, with details in App. C. This
approach is inspired by the unscreened perturbative bath expansion of Ref. 66,
but here adapted for a screened embedding.

Figure 3: Convergence of the CCSD local $2p_{z}$-orbital charge (top) and spin
(bottom) two-body fluctuations in $8\times 8$ Graphene sheets with increasing
size of embedded cluster, for both bare Coulomb and mRPA cluster interactions.
Inset: Error in these two-body correlators in $4\times 4$ $\textbf{k}-$point
meshes, where comparison to exact CCSD is possible.

We converge semi-metallic graphene sheets with a fully ab initio CCSD
description, embedding atomic fragments and systematically enlarging the bath
space for each fragment. On enlarging, the bath interactions provide increased
screening of the fragments explicitly, and the mRPA static screening is
correspondingly reduced in a fashion that precludes double-counting of the
fragment screening and naturally converges to an in-method (CCSD) exact limit.
Figure 3 shows fully ab initio two-body instantaneous charge and spin
fluctuations in the local atomic $2p_{z}$ space, on top of a symmetry-
preserving Hartree–Fock reference [69, 67]. The correlated treatment
quantitatively changes these two-body correlators, but bare subspace
interactions overestimate the magnitude of the changes from mean-field.
Importantly, the bias in these correlators with the bare interactions appears
unable to be compensated for with increasing bath, indicating the importance
of coupling to truly long-range $\mathbf{q}\rightarrow 0$ plasmonic modes for
an appropriate relaxation of these local properties, impractical in a local
embedding with bare interactions.

In contrast, the mRPA screened interactions fold in this coupling, resulting
in a rapid and stable convergence with increasing bath size, corroborated with
smaller $4\times 4$ $\mathbf{k}$-point meshes where comparison to the full
system CCSD result is possible and the same qualitative behaviour is observed.
While $\mathcal{U}_{\mathrm{mRPA}}$ is static, it nevertheless integrates over
all RPA diagrams, including plasmonic contributions required to appropriately
relax the subspace. This even impacts on the convergence of magnetic
correlators in the subspace, despite not being described in the
$\mathbf{\eta}^{(0)}$ description of RPA density fluctuations. The results
indicate that the plasmonic coupling therefore suppresses the tendency towards
formation of atomic magnetic moments in graphene. In App. D we also
demonstrate improvements in non-local magnetic fluctuations as well as
energetics in this system, and the insensitivity of these results to choice of
underlying reference mean-field theory.

In summary, the developed ‘mRPA’ efficiently integrates over all external RPA
diagrams to provide a manifestly static effective interaction for low energy
or local models, conserving all instantaneous RPA correlators in a chosen
subspace by construction. For multi-method workflows, we show this provides a
systematic and principled route for the inclusion of long-range and high-
energy screening in a local correlation framework. This also opens
opportunities in the push for fully improvable and non-empirical quantum
embedding, with reliable convergence in extended systems coupled to long-range
coherent quasiparticles.

###### Acknowledgements.

The authors gratefully acknowledge support from the European Union’s Horizon
2020 research and innovation programme under grant agreement No. 759063. We
are grateful to the UK Materials and Molecular Modelling Hub for computational
resources, which is partially funded by EPSRC (EP/P020194/1 and EP/T022213/1).
All calculations were performed and can be reproduced with the Vayesta quantum
embedding code 222Vayesta code available at
https://github.com/BoothGroup/Vayesta..

## References

  * Hubbard and Flowers [1964] J. Hubbard and B. H. Flowers, Proc. Roy. Soc. London A 277, 237 (1964).
  * Jørgensen [1975] F. Jørgensen, Mol. Phys. 29, 1137 (1975).
  * Shavitt and Redmon [2008] I. Shavitt and L. T. Redmon, J. Chem. Phys. 73, 5711 (2008).
  * Löwdin [2004] P. Löwdin, J. Chem. Phys. 19, 1396 (2004).
  * Ma  _et al._ [2021] H. Ma, N. Sheng, M. Govoni, and G. Galli, J. Chem. Theory Comput. 17, 2116 (2021).
  * Zheng  _et al._ [2018] H. Zheng, H. J. Changlani, K. T. Williams, B. Busemeyer, and L. K. Wagner, Front. Phys. 6 (2018).
  * Chang and Wagner [2023] Y. Chang and L. K. Wagner, “Learning emergent models from ab initio many-body calculations,” (2023), arXiv:2302.02899 [cond-mat.str-el] .
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Rev. Mod. Phys. 90, 025003 (2018).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * Sun and Chan [2016] Q. Sun and G. K.-L. Chan, Acc. Chem. Res. 49, 2705 (2016).
  * Ditte  _et al._ [2023] M. Ditte, M. Barborini, L. M. Sandonas, and A. Tkatchenko, “Molecules in environments: Towards systematic quantum embedding of electrons and drude oscillators,” (2023), arXiv:2303.16250 [physics.chem-ph] .
  * Nomura  _et al._ [2012] Y. Nomura, M. Kaltak, K. Nakamura, C. Taranto, S. Sakai, A. Toschi, R. Arita, K. Held, G. Kresse, and M. Imada, Phys. Rev. B 86, 085117 (2012).
  * Şaşıoğlu  _et al._ [2011] E. Şaşıoğlu, C. Friedrich, and S. Blügel, Phys. Rev. B 83, 121101 (2011).
  * Springer and Aryasetiawan [1998] M. Springer and F. Aryasetiawan, Phys. Rev. B 57, 4364 (1998).
  * Kinza and Honerkamp [2015] M. Kinza and C. Honerkamp, Phys. Rev. B 92, 045113 (2015).
  * Han  _et al._ [2021] X.-J. Han, P. Werner, and C. Honerkamp, Phys. Rev. B 103, 125130 (2021).
  * Aryasetiawan  _et al._ [2009] F. Aryasetiawan, J. M. Tomczak, T. Miyake, and R. Sakuma, Phys. Rev. Lett. 102, 176402 (2009).
  * Park  _et al._ [2020] J. W. Park, R. Al-Saadon, M. K. MacLeod, T. Shiozaki, and B. Vlaisavljevich, Chem. Rev. 120, 5878 (2020).
  * Katukuri  _et al._ [2020] V. Katukuri, N. Bogdanov, O. Weser, J. Brink, and A. Alavi, Phys. Rev. B 102 (2020).
  * Chen  _et al._ [2017] G. P. Chen, V. K. Voora, M. M. Agee, S. G. Balasubramani, and F. Furche, Annu. Rev. Phys. Chem. 68, 421 (2017).
  * Jørgensen and Sauer [2020] M. W. Jørgensen and S. P. A. Sauer, J. Chem. Phys. 152, 234101 (2020).
  * Bohm and Pines [1953] D. Bohm and D. Pines, Phys. Rev. 92, 609 (1953).
  * Gell-Mann and Brueckner [1957] M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).
  * Shishkin and Kresse [2007] M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).
  * van Loon  _et al._ [2021] E. G. C. P. van Loon, M. Rösner, M. I. Katsnelson, and T. O. Wehling, Phys. Rev. B 104, 045134 (2021).
  * Aryasetiawan  _et al._ [2004] F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Phys. Rev. B 70, 195104 (2004).
  * Miyake and Aryasetiawan [2008] T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008).
  * Shinaoka  _et al._ [2015] H. Shinaoka, M. Troyer, and P. Werner, Phys. Rev. B 91, 245156 (2015).
  * Vaugier  _et al._ [2012] L. Vaugier, H. Jiang, and S. Biermann, Phys. Rev. B 86, 165105 (2012).
  * Hausoel  _et al._ [2017] A. Hausoel, M. Karolak, E. Şaş$\iota$oğlu, A. Lichtenstein, K. Held, A. Katanin, A. Toschi, and G. Sangiovanni, Nat. Comms. 8, 16062 (2017).
  * Sakuma and Aryasetiawan [2013] R. Sakuma and F. Aryasetiawan, Phys. Rev. B 87, 165118 (2013).
  * van Loon  _et al._ [2018] E. G. C. P. van Loon, M. Rösner, G. Schönhoff, M. I. Katsnelson, and T. O. Wehling, npj Quantum Materials 3, 32 (2018).
  * Werner  _et al._ [2015] P. Werner, R. Sakuma, F. Nilsson, and F. Aryasetiawan, Phys. Rev. B 91, 125142 (2015).
  * Jang  _et al._ [2016] S. W. Jang, H. Sakakibara, H. Kino, T. Kotani, K. Kuroki, and M. J. Han, Sci. Rep. 6, 33397 (2016).
  * Chen  _et al._ [2015] J. Chen, A. J. Millis, and C. A. Marianetti, Phys. Rev. B 91, 241111 (2015).
  * Nomura  _et al._ [2019] Y. Nomura, M. Hirayama, T. Tadano, Y. Yoshimoto, K. Nakamura, and R. Arita, Phys. Rev. B 100, 205138 (2019).
  * Martins  _et al._ [2011] C. Martins, M. Aichhorn, L. Vaugier, and S. Biermann, Phys. Rev. Lett. 107, 266404 (2011).
  * Shih  _et al._ [2012] B.-C. Shih, Y. Zhang, W. Zhang, and P. Zhang, Phys. Rev. B 85, 045132 (2012).
  * Yoshimi  _et al._ [2022] K. Yoshimi, T. Misawa, T. Tsumuraya, and H. Seo, “Comprehensive $ab$ $initio$ investigation of the phase diagram of quasi-one-dimensional molecular solids,” (2022), arXiv:2210.13726 [cond-mat.str-el] .
  * Scott and Booth [2021] C. J. C. Scott and G. H. Booth, Phys. Rev. B 104, 245114 (2021).
  * Sheng  _et al._ [2022] N. Sheng, C. Vorwerk, M. Govoni, and G. Galli, J. Chem. Theory Comput. 18, 3512 (2022).
  * Honerkamp  _et al._ [2018] C. Honerkamp, H. Shinaoka, F. F. Assaad, and P. Werner, Phys. Rev. B 98, 235151 (2018).
  * Casula  _et al._ [2012] M. Casula, P. Werner, L. Vaugier, F. Aryasetiawan, T. Miyake, A. J. Millis, and S. Biermann, Phys. Rev. Lett. 109, 126408 (2012).
  * Miyake  _et al._ [2009] T. Miyake, F. Aryasetiawan, and M. Imada, Phys. Rev. B 80, 155134 (2009).
  * Aryasetiawan  _et al._ [2011] F. Aryasetiawan, T. Miyake, and R. Sakuma, in  _The LDA+DMFT approach to strongly correlated materials , 2011 Reihe Modeling and Simulation, Vol. 1_ , edited by E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein (Verlag des Forschungszentrum Jülich, 2011) Chap. 7.
  * Casida [1995] M. E. Casida, in  _Recent Advances In Density Functional Methods: (Part I)_ (World Scientific, 1995) pp. 155–192.
  * Scott  _et al._ [2023] C. J. C. Scott, O. J. Backhouse, and G. H. Booth, J. Chem. Phys. 158, 124102 (2023).
  * Toulouse  _et al._ [2009] J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Angyán, Phys. Rev. Lett. 102, 096404 (2009).
  * Toulouse  _et al._ [2010] J. Toulouse, W. Zhu, J. G. Angyán, and A. Savin, Phys. Rev. A 82, 032502 (2010).
  * Ángyán  _et al._ [2011] J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. 7, 3116 (2011).
  * Knizia and Chan [2012] G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
  * Sekaran  _et al._ [2023] S. Sekaran, O. Bindech, and E. Fromager, “A unified density-matrix functional construction of quantum baths in density matrix embedding theory beyond the mean-field approximation,” (2023), arXiv:2304.14729 [physics.chem-ph] .
  * Bickers and Scalapino [1989] N. Bickers and D. Scalapino, Ann. Phys. 193, 206 (1989).
  * van Aggelen  _et al._ [2013] H. van Aggelen, Y. Yang, and W. Yang, Phys. Rev. A 88, 030501 (2013).
  * Gukelberger  _et al._ [2015] J. Gukelberger, L. Huang, and P. Werner, Phys. Rev. B 91, 235114 (2015).
  * Springer  _et al._ [1998] M. Springer, F. Aryasetiawan, and K. Karlsson, Phys. Rev. Lett. 80, 2389 (1998).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Comp. Mol. Sci. 8, e1340 (2018).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020).
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Note [1] Coupled-cluster is defined in terms of the Fock matrix and fluctuation potential. The virtual-virtual block of the Fock matrix is shifted by $\Delta f_{ab}=\sum@\slimits@_{i\in\textrm{occ}}(ia|bi)-\mathcal{U}_{ia,ib}$ due to the screened $ph$ mRPA interactions, but other blocks are unchanged ensuring the reference-independence of mRPA.
  * Karton  _et al._ [2011] A. Karton, S. Daon, and J. M. Martin, Chem. Phys. Lett. 510, 165 (2011).
  * Knizia [2013] G. Knizia, J. Chem. Theory Comput. 9, 4834 (2013).
  * Knizia and Chan [2013] G. Knizia and G. K.-L. Chan, J. Chem. Theory Comput. 9, 1428 (2013).
  * Wouters  _et al._ [2016] S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, J. Chem. Theory Comput. 12, 2706 (2016).
  * Cui  _et al._ [2020] Z.-H. Cui, T. Zhu, and G. K.-L. Chan, J. Chem. Theory Comput. 16, 119 (2020).
  * Nusspickel and Booth [2022] M. Nusspickel and G. H. Booth, Phys. Rev. X 12, 011046 (2022).
  * Nusspickel  _et al._ [2023] M. Nusspickel, B. Ibrahim, and G. H. Booth, J. Chem. Theory Comput. 19, 2769 (2023).
  * Wehling  _et al._ [2011] T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).
  * Pokhilko and Zgid [2021] P. Pokhilko and D. Zgid, J. Chem. Phys. 155, 024101 (2021).
  * Note [2] Vayesta code available at https://github.com/BoothGroup/Vayesta.
  * Eshuis  _et al._ [2010] H. Eshuis, J. Yarkony, and F. Furche, J. Chem. Phys. 132, 234114 (2010).
  * Hale  _et al._ [2008] N. Hale, N. J. Higham, and L. N. Trefethen, SIAM J. Numer. Anal. 46, 2505 (2008).
  * Yeh  _et al._ [2022] C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, Phys. Rev. B 106, 235104 (2022).
  * Zhu and Chan [2021] T. Zhu and G. K.-L. Chan, J. Chem. Theory Comput. 17, 727 (2021).
  * Kosov [2017] D. Kosov, Chem. Phys. Lett. 690, 20 (2017).
  * McClain  _et al._ [2017] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, J. Chem. Theory Comput. 13, 1209 (2017).

##  Appendix A Efficient computation of RPA subspace moments

We consider an efficient approach for the construction of the subspace
projection of the RPA quantities $\mathbf{\eta}^{(0)}$ and
$\mathbf{\eta}^{(1)}$. These are defined by their relation in Eq. 5, as
required for the subsequent computation of the subspace screened interaction
$\mathcal{U}_{\mathrm{mRPA}}$ via Eq. 6. This is similar to derivations found
in Ref. 47 in the context of efficient $GW$ theory, which itself is based on
efficient RPA correlation energy computation of Ref. 71. Ultimately, since the
working equations of mRPA below are algorithmically similar to many RPA and
$GW$ methods, we expect the computational feasibility of efficient mRPA
implementations to be similar. Starting from Eq. 5, we can multiply both sides
from the right by $(\mathbf{A}+\mathbf{B})$ to give

$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})=[\mathbf{\eta}^{(0)}(\mathbf{A}+\mathbf{B})]^{2}.$
(7)

We then employ a standard low-rank decomposition of the two-electron repulsion
integrals (e.g. via density fitting or Cholesky decomposition) in the
particle-hole channel, as

$v_{ia,jb}\simeq\sum_{P}^{N_{\mathrm{aux}}}V_{ia,P}V_{jb,P}=\mathbf{V}\mathbf{V}^{T},$
(8)

where we use $P,Q,\dots$ to index elements of this auxiliary (RI) basis, whose
dimension $N_{\mathrm{aux}}$ scales $\mathcal{O}\left(N\right)$ with system
size, and $i,j,\dots$ ($a,b,\dots$) denote hole (particle) spin-orbital
energies from the mean-field reference description (note App. B details direct
construction of the spin-integrated mRPA interaction). Writing
$(\mathbf{A}\pm\mathbf{B})$ explicitly in the (direct) RPA form, as

$\displaystyle(\mathbf{A}-\mathbf{B})_{ia,jb}$
$\displaystyle=(\epsilon_{a}-\epsilon_{i})\delta_{ij}\delta_{ab}=\mathbf{\Delta},$
(9) $\displaystyle(\mathbf{A}+\mathbf{B})_{ia,jb}$
$\displaystyle=(\epsilon_{a}-\epsilon_{i})\delta_{ij}\delta_{ab}+2\sum_{P}V_{ia,P}V_{jb,P}$
$\displaystyle=\mathbf{\Delta}+2\mathbf{V}\mathbf{V}^{T},$ (10)

we can cast Eq. 7 as

$\displaystyle\mathbf{\eta}^{(0)}$
$\displaystyle=[(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})]^{1/2}(\mathbf{A}+\mathbf{B})^{-1}$
(11)
$\displaystyle=(\mathbf{\Delta}^{2}+2\mathbf{\Delta}\mathbf{VV}^{T})^{1/2}(\mathbf{\Delta}+2\mathbf{V}\mathbf{V}^{T})^{-1}.$
(12)

Using the Woodbury matrix identity, we can rewrite the first part of this term
as

$\displaystyle\mathbf{T}$
$\displaystyle=(\mathbf{\Delta}+2\mathbf{V}\mathbf{V}^{T})^{-1}$ (13)
$\displaystyle=\mathbf{\Delta}^{-1}-2\mathbf{\Delta}^{-1}\mathbf{V}(\mathbf{I}+2\mathbf{V}^{T}\mathbf{\Delta}^{-1}\mathbf{V})^{-1}\mathbf{V}^{T}\mathbf{\Delta}^{-1}.$
(14)

This now only requires the inversion of the diagonal matrix,
$\mathbf{\Delta}$, and a matrix of dimension $N_{\mathrm{aux}}$, with the
overall $ov\times N_{\mathrm{aux}}$ matrix able to be constructed in
$\mathcal{O}[N_{\mathrm{aux}}^{3}+N_{\mathrm{aux}}^{2}ov]$ time, where $o$ and
$v$ are the number of hole and particle states respectively. The second part
of the expression can be computed efficiently using the definition of the
matrix square-root as an integration in the complex plane [72],

$\mathbf{M}^{\frac{1}{2}}=\frac{1}{\pi}\int_{-\infty}^{\infty}\left(\mathbf{I}-z^{2}(\mathbf{M}+z^{2}\mathbf{I})^{-1}\right)dz.$
(15)

This results in

$\displaystyle\mathbf{\eta}^{(0)}$
$\displaystyle=(\mathbf{\Delta}^{2}+2\mathbf{\Delta
V}\mathbf{V}^{T})^{\frac{1}{2}}\mathbf{T}$ (16)
$\displaystyle=\frac{1}{\pi}\int_{-\infty}^{\infty}\left(\mathbf{I}-z^{2}(\mathbf{\Delta}^{2}+2\mathbf{\Delta
VV}^{T}+z^{2}\mathbf{I})^{-1}\right)dz\mathbf{T}.$ (17)

We can apply the Woodbury matrix identity a second time to this integrand, and
simplify the notation by introducing the intermediates,

$\displaystyle\mathbf{F}(z)$
$\displaystyle=(\mathbf{\Delta}^{2}+z^{2}\mathbf{I})^{-1}$ (18)
$\displaystyle\mathbf{Q}(z)$
$\displaystyle=2\mathbf{V}^{T}\mathbf{F}(z)\mathbf{\Delta V},$ (19)

where $\mathbf{F}(z)$ is a diagonal matrix in the $ov$ space, and
$\mathbf{Q}(z)$ is a $N_{\mathrm{aux}}\times N_{\mathrm{aux}}$ matrix which
can be constructed in $\mathcal{O}[N_{\mathrm{aux}}^{2}ov]$ time (the leading-
order computational step). This casts Eq. 17 into the form

$\displaystyle\eta^{(0)}=\frac{1}{\pi}$
$\displaystyle\int_{-\infty}^{\infty}\left[\mathbf{I}-z^{2}\mathbf{F}(z)+\right.$
(20) $\displaystyle\left.2z^{2}\mathbf{F}(z)\mathbf{\Delta
V}(\mathbf{I}+\mathbf{Q}(z))^{-1}\mathbf{V}^{T}\mathbf{F}(z)\right]dz\mathbf{T}.$

This integrand bears much similarity to quantities computed in efficient RPA
implementations, or Matsubara $GW$ implementations [73, 74, 71], where $z$
represents an imaginary frequency variable. Following Ref. 47, it is possible
to analytically integrate the mean-field contribution to this quantity, as
well as partially integrate the slowest order decaying contribution to the
integrand in $z$, to improve the efficiency of numerical quadrature. This
results in the expression

$\displaystyle\eta^{(0)}$ $\displaystyle=\mathbf{\Delta}\mathbf{T}$
$\displaystyle+2\int_{0}^{\infty}e^{-t\mathbf{\Delta}}\mathbf{\Delta
V}\mathbf{V}^{T}e^{-t\mathbf{\Delta}}\mathbf{T}dt$
$\displaystyle+\frac{2}{\pi}\int_{-\infty}^{\infty}z^{2}\mathbf{F}(z)\mathbf{\Delta
V}\left((\mathbf{I}+\mathbf{Q}(z))^{-1}-\mathbf{I}\right)\mathbf{V}^{T}\mathbf{F}(z)\mathbf{T}dz.$
(21)

The first numerical integral in Eq. 21 (where the integrand decays
exponentially) is computed via Gauss-Laguerre quadrature, while the second
(where the integrand decays as $\mathcal{O}[z^{-4}]$) is evaluated via
Clenshaw–Curtis quadrature. Quadrature points and weights are optimized
according to the protocol in Ref. 47, resulting in exponential convergence
with number of points, and a practical necessity for no more than $24$ points.

While this describes the calculation of the full system $\mathbf{\eta}^{(0)}$,
only the component of $\mathbf{\eta}^{(0)}$ in the subspace of excitations is
required. The projector from the space of full system particle-hole
excitations into this subspace is constructed as

$Z_{ia,pq}=C_{pi}C_{qa}=\mathbf{Z},$ (22)

where $C_{pi}$ are the coefficients of the model space in the basis of the
hole and particle reference states. Including $\mathbf{Z}$ explicitly on both
sides of Eq. 21 and optimizing the contraction order of all quantities, it can
be shown that the scaling for the construction of the subspace component is
$\mathcal{O}\left(N^{4}\right)$, with the rate limiting step the
$\mathcal{O}\left(N_{\textrm{aux}}ov\right)$ construction of the
$\mathbf{Q}(z)$ matrix. This can be reduced to $\mathcal{O}\left(N^{3}\right)$
with the introduction of doubly-factorized integrals such as those from tensor
hypercontraction or forms found in the space-time approach to $GW$ theory, or
from an iterative approach to the direct construction and inversion of the
$(\mathbf{I}+\mathbf{Q}(z))$ matrix with local projectors, which will be
considered in future work. Furthermore, since the $\mathbf{\eta}^{(1)}$
quantity is diagonal, finding the subspace component of this quantity is
trivial.

Finally, we note that this $\mathcal{O}\left(N^{4}\right)$ formal scaling is
unchanged if you want to find the mRPA screened interactions of a single
subspace (even one whose size grows linearly with system size) or multiple
local subspaces simultaneously, since the full system $\mathbf{Q}(z)$ matrix
construction and inversion is still performed once, and other contractions do
not scale worse than this. This feature is exploited in the main text in the
results of Fig. 2, where each atom is considered part of its own screened
subspace interactions, and thus the number of subspaces which you need to find
the screened interactions for increases linearly with system size.
Nevertheless, the formal scaling of finding all mRPA screened subspace
interactions is still $\mathcal{O}\left(N^{4}\right)$ overall.

##  Appendix B Spin independence of mRPA interactions

In the case of a spin-conserving reference state where the irreducible RPA
polarizability is the same for both spin-channels, and where the interaction
kernel is spin-independent (e.g. in the absence of magnetic fields), the mRPA
interaction resulting from the above procedure is also spin-independent. This
allows for straightforward compatibility with methods to solve the effective
Hamiltonian designed for use with the bare Coulomb interaction in ‘restricted’
formulations, and furthermore allows for efficiency improvements in the
construction of this effective interaction.

In such a case, $\mathbf{V}\mathbf{V}^{T}$ is the same for all spin-conserving
transitions, so the Coulomb coupling between two spin-conserving excitations
depends only on the spatial orbitals. This combined with the fact that
$\mathbf{\Delta}$ is spin-blocked means that when combining Eqs. 14 and 21 all
terms containing the expansion of the Coulomb interaction over an auxiliary
space intermediate are spin-independent. The only spin-dependent contribution
is thus the product of the two non-Coulombic terms in these equations,
$\mathbf{\Delta}\mathbf{\Delta}^{-1}=\mathbf{I}$, and all other terms must be
spin-independent.

When separated into blocks of single-particle excitations in each
$(\alpha,\beta)$ spin channel the dRPA zeroth dd-moment must as such take the
form

$\mathbf{\eta}^{(0)}=\left(\begin{aligned} &\mathbf{I}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{I}\\\ \end{aligned}\right)+\left(\begin{aligned}
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
\end{aligned}\right),$ (23)

where $\mathbf{\eta}^{(0)}_{\text{s}}$ is a spin-independent, correlation-
induced part of the zeroth dd-moment, depending only on the spatial component
of the orbital basis. In this case, we expand the inverse of
$\mathbf{\eta}^{(0)}$ in the spin-orbital basis via a Neumann series, as

$\displaystyle(\mathbf{\eta}^{(0)})^{-1}$
$\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}\left(\begin{aligned}
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
\end{aligned}\right)^{n}$ (24) $\displaystyle=\left(\begin{aligned}
&\mathbf{I}&\mathbf{0}\\\ &\mathbf{0}&\mathbf{I}\\\
\end{aligned}\right)-\left(\begin{aligned}
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
\end{aligned}\right)\sum_{n=1}^{\infty}(-2)^{n-1}\left(\begin{aligned}
&(\mathbf{\eta}^{(0)}_{\text{s}})^{n-1}&(\mathbf{\eta}^{(0)}_{\text{s}})^{n-1}\\\
&(\mathbf{\eta}^{(0)}_{\text{s}})^{n-1}&(\mathbf{\eta}^{(0)}_{\text{s}})^{n-1}\\\
\end{aligned}\right)$ (25) $\displaystyle=\left(\begin{aligned}
&\mathbf{I}&\mathbf{0}\\\ &\mathbf{0}&\mathbf{I}\\\
\end{aligned}\right)-\left(\begin{aligned}
&\mathbf{\eta}^{(0)}_{\text{s}}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{\eta}^{(0)}_{\text{s}}\\\
\end{aligned}\right)\left(\begin{aligned}
&(\mathbf{I}+2\mathbf{\eta}^{(0)}_{\text{s}})^{-1}&(\mathbf{I}+2\mathbf{\eta}^{(0)}_{\text{s}})^{-1}\\\
&(\mathbf{I}+2\mathbf{\eta}^{(0)}_{\text{s}})^{-1}&(\mathbf{I}+2\mathbf{\eta}^{(0)}_{\text{s}})^{-1}\\\
\end{aligned}\right)$ (26)

We can therefore write the inverse of the zeroth dd-moment compactly as

${\mathbf{\eta}^{(0)}}^{-1}=\left(\begin{aligned} &\mathbf{I}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{I}\\\ \end{aligned}\right)-\left(\begin{aligned}
&\mathbf{\kappa}&\mathbf{\kappa}\\\ &\mathbf{\kappa}&\mathbf{\kappa}\\\
\end{aligned}\right)$ (27)

with

$\mathbf{\kappa}=\mathbf{\eta}^{(0)}_{\text{s}}\left(1+2\mathbf{\eta}^{(0)}_{\text{s}}\right)^{-1}.$
(28)

Insertion of this definition into Eq. 6 gives the spin-independent mRPA
interactions for restricted references as

$v_{\text{s}}=\mathbf{\kappa}\mathbf{\Delta}\mathbf{\kappa}-\frac{1}{2}\left(\mathbf{\kappa}\mathbf{\Delta}+\mathbf{\Delta}\mathbf{\kappa}\right).$
(29)

Given that this only requires an inversion of a matrix in Eq. 28 which is half
the size of the analogous spinned quantity in Eq. 6, this can result in a more
computationally efficient implementation with respect to model space size,
where applicable.

##  Appendix C RPA interacting bath orbitals

We seek to define a bath expansion for a quantum embedding to systematically
expand the single particle space, in an optimally efficient way to describe
the physics of the zeroth moment of the dd-response at the RPA level
($\mathbf{\eta}^{(0)}_{ia,jb}$). From a viewpoint of the static expectation
values, this can be considered as the optimal single-particle bath expansion
to span the fragment correlations of the two-body reduced density matrix at
the RPA level. While the use of mRPA interactions in the cluster ensures that
the subspace $\mathbf{\eta}^{(0)}$ is exactly reproduced at the level of RPA,
this is not true when the cluster is solved at a higher level of theory, and
thus a bath expansion can provide a route to systematically control the
remaining error in this beyond-RPA description of the cluster.

We define a bath expansion which will optimally minimize the error in the
cluster $\mathbf{\eta}^{(0)}$ with bare interactions, ensuring that the bath
spans these most relevant degrees of freedom at the level of the resulting
solver, and in the process minimize the difference between the mRPA and bare
interactions throughout the cluster for a given number of bath orbitals. This
approach is related to a similar (interacting) bath expansion which was
introduced in Ref. [66] which optimally captured the one-body physics from a
bare (Møller-Plesset) second-order perturbative level of theory. We expect
that a screened perturbative bath expansion based on RPA to be more suitable
for extended and polarizable systems, while noting that its construction
relies on knowledge of the zeroth moment of the dd-response which is already
required to construct the mRPA interactions in the resulting cluster, making
the approach particularly straightforward and consistent in this framework.

We start from a cluster formed from an arbitrary set of (generally local)
‘fragment’ orbitals. These are augmented with the bath space from density
matrix embedding theory (DMET), which ensures that the reference single
determinant state is exactly spanned by a single state in the cluster, and
that the mean-field density matrix of the fragment can be faithfully
reproduced [51, 63, 52, 64, 65, 66]. We note that this construction also
ensures the applicability of mRPA directly to this cluster. To obtain the RPA
bath expansion orbitals to further augment this initial cluster, we form one-
particle/hole pseudo-density matrices representing the correlated hole and
particle fluctuations at the level of RPA which couple this initial cluster to
its wider environment,

$\displaystyle\gamma^{\mathbf{x}}_{ij}$
$\displaystyle=-\sum\limits_{\tilde{a}}(\mathbf{\eta}^{(0)}_{i\tilde{a}j\tilde{a}}-\delta_{ij}\delta_{\tilde{a}\tilde{a}})$
(30) $\displaystyle\gamma^{\mathbf{x}}_{ab}$
$\displaystyle=-\sum\limits_{\tilde{i}}(\mathbf{\eta}^{(0)}_{\tilde{i}a\tilde{i}b}-\delta_{ab}\delta_{\tilde{i}\tilde{i}}).$
(31)

Here, $\tilde{i}$, ($\tilde{a}$) indices denote the local hole (particle)
states associated with the initial cluster, $\mathbf{x}$, which can be
efficiently implemented via projection operators, while other indices
represent environmental hole (particle) orbitals in the complementary space
[66]. Note that all indices represent spin-orbitals, which in spin colinear
systems will require that all indices have the same spin label. As such, from
this point on all indices will be assumed to represent the same spin channel.

The eigenvectors of these matrices represent orbitals of a bath expansion of
the correlated environment in occupied and unoccupied states, with their
corresponding eigenvalues quantifying their importance in describing the full
RPA physics of $\mathbf{\eta}^{(0)}$ in the cluster. These eigenvalues can
then be used as a systematic truncation via a threshold in order to define the
size of the bath space, analogous to the approach in Ref. 66. The mRPA (or
bare) interactions can then be found in this subspace for the high-level
solver in a quantum embedding framework. The required
$\mathbf{\eta}^{(0)}_{i\tilde{a}j\tilde{a}}$ components can be calculated in
$\mathcal{O}\left(N^{4}\right)$ computational cost based on the approach in
App. A, with additional techniques (not pursued here) expected to be able to
reduce this to $\mathcal{O}\left(N^{3}\right)$ via restriction of the full
excitation space in the construction of
$\mathbf{\eta}^{(0)}_{i\tilde{a}j\tilde{a}}$, or working with a doubly
factorized bare interaction [47].

To provide some further insight into the nature of these bath orbitals, we can
explicitly consider the relation between the two-body density matrix,
$\Gamma_{pqrs}$, and the zeroth moment of the four-point dd-response in an
arbitrary basis as an expansion over all neutral excitations,

$\chi^{(0)}_{pqrs}=\sum_{\nu\neq
0}\bra{0}\hat{c}_{q}^{\dagger}\hat{c}_{p}\ket{\nu}\bra{\nu}\hat{c}_{s}^{\dagger}\hat{c}_{r}\ket{0}=\Gamma_{pqsr}-\gamma_{pq}\gamma_{rs}+\delta_{qs}\gamma_{pr},$
(32)

where $\chi^{(0)}_{pqrs}$ is the zeroth moment of the dd-response, which is
related to the $\mathbf{\eta}^{(0)}$ quantity central to this work via
summation over the excitations and de-excitation manifold,

$\mathbf{\eta}^{(0)}_{iajb}=\chi^{(0)}_{iajb}+\chi^{(0)}_{iabj}+\chi^{(0)}_{aijb}+\chi^{(0)}_{aibj}.$
(33)

We find that the RPA pseudo-density matrices for each cluster $\mathbf{x}$
used to define the RPA bath expansion in Eqs. 30-31 are related to a
hypothetical RPA one-body density matrix, $\gamma$, as

$\displaystyle\gamma^{\mathbf{x}}_{ij}$
$\displaystyle=\sum_{\tilde{a}}\left[(\delta_{ij}-\gamma_{ij})(\delta_{\tilde{a}\tilde{a}}-\gamma_{\tilde{a}\tilde{a}})+\gamma_{\tilde{a}\tilde{a}}\gamma_{ij}+4\gamma_{i\tilde{a}}\gamma_{j\tilde{a}}\right]$
(34) $\displaystyle\gamma^{\mathbf{x}}_{ab}$
$\displaystyle=\sum_{\tilde{i}}\left[(\delta_{\tilde{i}\tilde{i}}-\gamma_{\tilde{i}\tilde{i}})(\delta_{ab}-\gamma_{ab})+\gamma_{\tilde{i}\tilde{i}}\gamma_{ab}+4\gamma_{\tilde{i}a}\gamma_{\tilde{i}b}\right].$
(35)

While this does not guarantee positive definiteness of the quantities in Eqs.
30-31, since the RPA density matrices are not well-defined [75] (hence the use
of pseudo-density matrix), it gives some physical justification for such a
property. It also demonstrates that the expansion captures the one-body
coupling to the environment, despite originating from a purely two-body
correlated quantity.

##  Appendix D Non-local and energetic expectation values in graphene

Figure 4: Analogous results to Fig. 3 of the main text, but all bath orbitals
and mRPA interactions are constructed on a Kohn–Sham PBE rather than
Hartree–Fock reference state. Convergence of the CCSD local $2p_{z}$-orbital
charge (top) and spin (bottom) two-body fluctuations in $8\times 8$ Graphene
with increasing size of embedded cluster, for both bare Coulomb and mRPA
cluster interactions. Inset: Error in these two-body correlators in $4\times
4$ $\textbf{k}-$point meshes, where comparison to exact CCSD (on the PBE
reference) is possible.

In this section, we expand on the results of mRPA and quantum embedding
calculations on the graphene system, considering the effect of different
reference states, as well as the impact of the mRPA screening on non-local
expectation values. The results of the main text were performed for a $4\times
4$ and $8\times 8$ $\Gamma$-centered $\mathbf{k}$-point mesh, in a cc-pVDZ
basis with a restricted Hartree–Fock reference state, via the Vayesta quantum
embedding package, which interfaces with the pyscf codebase [58, 57, 65, 76].
However, no double-counting results from the use of Kohn–Sham references in
order to define the orbitals and their energies in the irreducible
polarizability instead of Hartree–Fock. There may be advantages in this in
terms of convergence of the mean-field solution in metallic systems and
accuracy of the resulting mRPA interactions in extended systems. However, too
much dependence on this starting reference would also be a cause of concern
and an indication of also inheriting the empiricism in the choice of mean-
field theory.

In Fig. 4 we present analogous calculations to those of Fig. 3 in the main
text, showing the local charge and spin correlators with bare and mRPA
interactions in the quantum embedding of increasing cluster sizes, but now on
a Kohn–Sham reference state with the PBE functional. This can change both the
construction of the RPA bath orbitals of App. C defining the single-particle
space of the cluster, as well as the mRPA interactions which are used in its
subsequent solution. While the mean-field gap is quite different to the
Hartree–Fock reference, the convergence of the CCSD local correlators from the
clusters for both mRPA and bare interactions is almost identical. This is
reassuring that it is genuinely describing the effect of the post-mean-field
physics, and inherits the largely reference-independent property of the CCSD
solver in the full system limit. Future work will consider a self-consistent
approach to the mRPA interactions, formally removing the mean-field
dependence.

Figure 5: Convergence of the CCSD nearest-neighbor $2p_{z}$-orbital spin
correlated fluctuations in $8\times 8$ Graphene with increasing size of
embedded cluster, for both bare Coulomb and mRPA cluster interactions. Inset:
Error in this two-body correlators in $4\times 4$ $\textbf{k}-$point mesh,
where comparison to exact CCSD is possible.

We also consider the accuracy of non-local expectation values in a quantum
embedding framework with mRPA screened interactions, where the operators span
more than one fragment (impurity). One approach is to simply enlarge the
fragment space to encompass the range of the operator of interest. However,
this is inefficient and will not be suitable for full system expectation
values such as total energies. Reference 67 considered approaches to couple
the solutions of different embedded quantum clusters for non-local expectation
values. We therefore follow this prescription to compute atomic $2p_{z}$ spin
correlators to quantify the fluctuations between two nearest-neighbour atoms
(labelled $A$ and $B$), as $\langle{\hat{S}_{z}^{A}}{\hat{S}_{z}^{B}\rangle}$
in Fig. 5, via construction of the two-body reduced density matrix (denoted
$\Gamma[(\gamma,K^{*})[\Psi^{\mathbf{x}}]]$ in Ref. 67). We highlight that the
atomic fragmentation means that these correlators span different embedded
cluster solutions.

These results again show the impact of the mRPA screening in the cluster
solutions for this system, even for non-local expectation values which are not
explicitly conserved in the mRPA construction. Plasmonic modes over the
longest length-scales are folded in to the local interactions, reducing the
nearest-neighbour magnetic fluctuations, and enabling high-accuracy and
systematic convergence compared to bare interacting clusters.

When extending this analysis to total energies, some more care is required.
The magnetic correlators considered extend beyond the length of any single
fragment, but the union of all fragments do span the considered operators.
This is not always the case however, and in our embedding protocol where
fragments are intrinsic atomic orbitals of minimal size, the space of high-
energy unoccupied states in large basis sets is not contained in any fragment.
While the mRPA interactions will relax and improve the description of the
correlated low-energy fragment spaces when solved in isolation, they will not
account for the incomplete trace in the evaluation of non-local expectation
values over these high energy degrees of freedom not spanned by any of the
solved subspaces.

When considering total energies with variational solvers of fragmented (yet
incomplete) subspaces, this incompleteness in the space spanned by the
Hamiltonian in the union of the subspaces will serve to increase the total
energy of the fragmented system. Similarly, the use of mRPA interactions in
the subspaces will serve to also increase the energy of the individual
subspace descriptions. Despite the fact that the description of these
subspaces with mRPA interactions is closer to the projection of the full
system description into this subspace, the use of mRPA interactions will
therefore necessarily serve to increase the total energy error. This is shown
in Fig. 6, where the correlation energy error of graphene increases with mRPA
interactions. This energy over the embedded fragments is defined via the
‘projected’ energy estimate, introduced in Ref. 66, and expanded on in Ref.
67. It should be noted that this argument of energy error increase with mRPA
interactions is not necessarily true when considering energy differences.

Figure 6: Convergence of the fragmented CCSD correlation energy in $8\times 8$
Graphene with increasing size of embedded cluster, for bare Coulomb
interactions and mRPA cluster interactions with and without the $\Delta E_{k}$
RPA energy correction. Inset: Equivalent error in the correlation energy in a
$4\times 4$ $\textbf{k}-$point mesh, where comparison to full system CCSD is
possible. Calculations were based upon a RHF reference.

We therefore seek a way to improve the convergence of energy estimators in
keeping with the spirit of mRPA, to include the energetic contributions from
high-energy and long-range contributions that are not contained within any
subspace description. This would therefore directly include the energetic
effect of e.g. plasmon modes, rather than just the effect of their coupling on
the low-energy descriptions. The brute-force approach via an (interacting)
bath expansion can be seen in Fig. 6 to be poorly convergent, while previous
DMET approaches have considered direct enlargement of fragment spaces to
include these energetic contributions, but substantially increases the cost of
the subspace solutions [65].

As an alternative approach, we formulate a full-system energetic correction to
these incomplete fragmented subspace expectation values via the mRPA
construction. We begin by noting that we can construct a global two-body
density matrix contribution from the RPA via the relation given in Eq. 32.
$\chi^{(0)}_{pqrs}$ then relates to the cumulant (or connected) component of
the two-body density matrix, $K_{pqrs}$, as

$K_{pqrs}=\chi^{(0)}_{pqsr}-\gamma_{ps}(\delta_{rq}-\gamma_{rq}).$ (36)

This represents the correlated, two-body correlations not decomposable as a
combination of one-body correlations. The contribution of these irreducible
two-body fluctuations to the total energy is given by

$\displaystyle E_{k}$ $\displaystyle=\frac{1}{2}\sum_{pqrs}(pq|rs)K_{pqrs}$
(37)
$\displaystyle=\frac{1}{2}\sum_{pqrs}(pq|rs)\left(\chi^{(0)}_{pqsr}-\gamma_{ps}(\delta_{rq}-\gamma_{rq})\right).$
(38)

Under the RPA description of $ph$ (de-)excitations, the first term naturally
becomes

$\displaystyle\sum_{pqrs}(pq|rs)\chi^{(0)}_{pqsr}=\sum_{\begin{subarray}{c}ij\in\text{occ}\\\
ab\in\text{vir}\end{subarray}}(ia|jb)\mathbf{\eta}^{(0)}_{iajb}$ (39)

while approximating $\gamma$ by its mean-field value in the remaining terms
gives

$\displaystyle\sum_{pqrs}(pq|rs)\gamma_{ps}(\delta_{rq}-\gamma_{rq})$
$\displaystyle\approx\sum_{pqrs}\delta_{ps}^{p\in\text{occ}}\delta_{rq}^{r\in\text{vir}}(pq|rs)$
$\displaystyle=\sum_{\begin{subarray}{c}i\in\text{occ}\\\
a\in\text{vir}\end{subarray}}(ia|ia).$ (40)

We therefore obtain an overall correlation energy contribution from these
irreducible two-body fluctuations as

$E_{k}=\frac{1}{2}\text{Tr}\left[\left(\mathbf{\eta}^{(0)}-\mathbf{I}\right)\mathbf{v}\right].$
(41)

To use this as a correction to our correlated energy estimate obtained via
quantum embedding of our incomplete fragmented wavefunction, we must remove
local contributions from each cluster to avoid double counting. This can be
achieved simply in our quantum embedding approach by computing a contribution
to these two-body energetic fluctuations in each subspace, by projecting one
index into the fragment space, and the remaining indices restricted to the
subspace of consideration [67]. This provides a decomposition of the total RPA
two-body energy for a cluster $\mathbf{x}$ as

$E_{k}^{\mathbf{x}}=\frac{1}{2}\text{Tr}\left[\mathbf{P}^{\mathbf{x}}\left(\mathbf{\eta}^{(0)}-\mathbf{I}\right)\mathbf{v}\right],$
(42)

where $\mathbf{P}^{\mathbf{x}}$ is defined as symmetrically projecting one of
the occupied indices of the full-system $\mathbf{\eta}^{(0)}$ into the
fragment space, with all other indices projected into the cluster space of
$\mathbf{x}$ (as in Eq. 25 of [67]).

The final non-local energy correction is therefore obtained by deducting the
local subspace contributions from the global expectation, giving the result

$\Delta E_{k}=E_{k}-\sum_{\mathbf{x}}E_{k}^{\mathbf{x}}.$ (43)

This is equivalent to only including contributions in Eq. 41 where all indices
are not contained within any of their cluster space. While this choice of
decomposition into subspace contributions to avoid double-counting is specific
for our approach for computing expectation values, the general approach could
be easily adapted for other embedding or subspace methodologies in order to
define a non-local RPA energy correction.

Application of this RPA energetic two-body correction for the incompleteness
of the cluster spaces in the atomic IAO quantum embedding of graphene is shown
in Fig. 6. This is found to substantially improve the rate of convergence with
respect to the size of the explicit CCSD clusters towards the full CCSD limit,
with the energy correction accurately compensating for the long-range bubble
contributions not spanned by the clusters. Other energy contributions are
expected to decay more rapidly, and therefore more easily described by the
increasing local cluster spaces.


# Eigenvalue Bounds for Saddle-Point Systems with Singular Leading Blocks

Susanne Bradley Department of Computer Science, The University of British
Columbia, Vancouver, Canada V6T 1Z4<EMAIL_ADDRESS>greif@cs.ubc.ca). Chen
Greif22footnotemark: 2

###### Abstract

We derive bounds on the eigenvalues of saddle-point matrices with singular
leading blocks. The technique of proof is based on augmentation. Our bounds
depend on the principal angles between the ranges or kernels of the matrix
blocks. Numerical experiments validate our analytical findings.

##  1 Introduction

Consider the saddle-point system

$\begin{bmatrix}A&B^{T}\\\ B&0\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}f\\\ g\end{bmatrix},$ (1)

where $A\in\mathbb{R}^{n\times n}$ is symmetric positive semidefinite and
$B\in\mathbb{R}^{m\times n}$ has full row rank, with $m<n$. We denote the
coefficient matrix by

$\mathcal{K}=\begin{bmatrix}A&B^{T}\\\ B&0\end{bmatrix}.$ (2)

We assume throughout that $\mathcal{K}$ is invertible. Our goal in this paper
is to derive eigenvalue bounds for $\mathcal{K}$ under the assumption that $A$
is singular.

#### Contribution of this paper

We derive a nonzero bound on the positive eigenvalues of $\mathcal{K}$ that
does not require invertibility of $A$, by considering the principal angles
between the ranges/kernels of $A$ and $B$. Ruiz et al. [11] also developed a
lower positive eigenvalue bound using principal angles, but their analysis
assumes a positive definite $A$.

#### Notation

Our analysis will rely on the eigenvalues and singular values of $A$ and $B$,
as well as some other matrices we will introduce later in the text. We will
denote the eigenvalues of a matrix $M\in\mathbb{R}^{n\times n}$ by

$\mu_{i}(M),\quad i=1,\dots,n,$

and in terms of ordering we will assume that

$\mu_{1}(M)\geq\mu_{2}(M)\geq\cdots\geq\mu_{n}(M).$

We follow the same convention for singular values of a rectangular matrix $N$,
but we use $\sigma$ rather than $\mu$: i.e., the the singular values of
$N\in\mathbb{R}^{m\times n}$ are denoted by

$\sigma_{1}(N)\geq\sigma_{2}(N)\geq\cdots\geq\mu_{m}(N)\geq 0.$

To increase clarity, we will often refer to the maximal eigenvalues/singular
values $\mu_{1}(M)$ and $\sigma_{1}(N)$ by $\mu_{\max}(M)$ and
$\sigma_{\max}(M)$ respectively. Similarly, we will refer to the minimal
values $\mu_{n}(M)$, $\sigma_{m}(N)$ by $\mu_{\min}(M)$ and
$\sigma_{\min}(M)$. The positive eigenvalues of a matrix will be denoted by a
“+” superscript – for instance, we denote the smallest nonzero eigenvalue of a
semidefinite matrix $M$ by $\mu_{\min}^{+}(M)$. For simplicity, we will omit
the arguments to $\mu$ and $\sigma$ when we refer to the eigenvalues and
singular values of $A$ and $B$. That is, we let

$\displaystyle\mu_{\max}$ $\displaystyle=\mu_{\max}(A);\ \ \sigma_{\max}$
$\displaystyle=\sigma_{\max}(B);$ $\displaystyle\mu_{\min}$
$\displaystyle=\mu_{\min}(A);\ \ \sigma_{\min}$
$\displaystyle=\sigma_{\min}(B);$ $\displaystyle\ \ \ \ \
\mu_{\min}^{+}=\mu_{\min}^{+}(A).$

#### Outline

In Section 2 we discuss our general approach of augmenting the leading block
of a saddle-point matrix to obtain a lower bound on the positive eigenvalues.
In Section 3 we provide new bounds, which rely on the angles between the
kernel of $A$ and $B$. We then present numerical experiments in Section 4 and
concluding remarks in Section 5.

##  2 Lower positive eigenvalue bounds using leading block augmentation

To illustrate the challenge posed by the problem in hand, recall the following
result of Rusten and Winther [12, Lemma 2.1]. In their analysis it is assumed
that $A$ is positive definite (as opposed to semidefinite); however, the proof
of this lemma does not rely on this, so the result still holds when $A$ is
semidefinite.

######  Lemma 1.

Then, the eigenvalues of $\mathcal{K}$ are bounded in the union of intervals

$I^{-}\ \cup\ I^{+},$

where

$I^{-}=\left[\frac{1}{2}(\mu_{\min}-\sqrt{\mu_{\min}^{2}+4\sigma_{max}^{2}}),\frac{1}{2}(\mu_{\max}-\sqrt{\mu_{\max}^{2}+4\sigma_{\min}^{2}})\right]$

and

$I^{+}=\left[\mu_{\min},\frac{1}{2}(\mu_{\max}+\sqrt{\mu_{\max}^{2}+4\sigma_{\max}^{2}})\right].$

When $A$ is singular, the upper bounds on both positive and negative values of
$\mathcal{K}$ are unchanged, and the lower negative bound reduces to
$-\sigma_{\max}$. The main difficulty is that the lower bound on positive
eigenvalues reduces to zero, which is not a useful bound, especially in
situations where $\mathcal{K}$ is known to be nonsingular (which is our
assumption throughout this paper). When the null spaces of $A$ and $B$ are
well separated, the matrix $\mathcal{K}$ may in fact be well-conditioned and
its minimal positive eigenvalue bounded away from zero.

As a motivating example that illustrates the range of possibilities, consider
the coefficient matrix

$\mathcal{K}=\begin{bmatrix}1&0&b_{1}\\\ 0&0&b_{2}\\\
b_{1}&b_{2}&0\end{bmatrix}\textrm{ where }A=\begin{bmatrix}1&0\\\
0&0\end{bmatrix}\textrm{ and }B=\begin{bmatrix}b_{1}&b_{2}\end{bmatrix},$ (3)

with $b_{1}^{2}+b_{2}^{2}=1$ and $b_{1},b_{2}>0$. The eigenvalues of $A$ and
singular value of $B$ are the same for all such $b_{1}$, $b_{2}$, but the
lowest positive eigenvalue of $\mathcal{K}$ varies depending on $b_{1}$ and
$b_{2}$. The eigenvalues $\lambda$ of $\mathcal{K}$ are the roots of the cubic
polynomial $p(\lambda)=\lambda^{3}-\lambda^{2}-\lambda+b_{2}^{2}$. This
polynomial has two positive roots and and one negative root [2, Corollary
2.2]; the smaller positive root approaches zero as $b_{2}$ goes to zero (i.e.,
when $A$ and $B$ have overlapping null spaces), but as $b_{2}$ goes to 1
(i.e., when $A$ and $B$ have orthogonal null spaces) the smaller positive root
approaches 1.

We now present a general approach for deriving nonzero bounds for the lower
positive eigenvalues of $\mathcal{K}$ when $A$ is singular. We recall the
following result [5, 6]:

######  Lemma 2.

Let

$\mathcal{K}(W)=\begin{bmatrix}A+B^{T}WB&B^{T}\\\ B&0\end{bmatrix},$ (4)

where $W\in\mathbb{R}^{m\times m}$. If $\mathcal{K}$ and $\mathcal{K}(W)$ are
both nonsingular, then

$\mathcal{K}^{-1}=(\mathcal{K}(W))^{-1}+\begin{bmatrix}0&0\\\
0&W\end{bmatrix}.$ (5)

We will assume that $W$ is positive semidefinite and the leading block
$A_{W}:=A+B^{T}WB$ of $\mathcal{K}(W)$ is positive definite. We can use this
along with (5) to derive a nonzero bound on the lower positive eignvalues of
$\mathcal{K}$, using a free matrix parameter $W$.

######  Theorem 3.

Let $W\in\mathbb{R}^{m\times m}$ be a symmetric positive semidefinite matrix
and let $A_{W}=A+B^{T}WB$. Then the positive eigenvalues of $\mathcal{K}$ are
greater than or equal to

$\min\left\\{\mu_{\min}(A_{W}),\frac{1}{\mu_{\max}(W)}\right\\}.$

###### Proof.

We derive a lower bound on the positive eigenvalues of $\mathcal{K}$ by
considering an upper bound on the eigenvalues of $\mathcal{K}^{-1}$. By
combining [1, Equation (3.4)] and (5), we obtain

$\mathcal{K}^{-1}=\begin{bmatrix}A_{W}^{-1}-A_{W}^{-1}B^{T}S_{W}^{-1}BA_{W}^{-1}&A_{W}^{-1}B^{T}S_{W}^{-1}\\\
S_{W}^{-1}BA_{W}^{-1}&-S_{W}^{-1}+W\end{bmatrix},$ (6)

where $S_{W}=BA_{W}^{-1}B^{T}$. Notice that we can write

$\mathcal{K}^{-1}=\begin{bmatrix}A_{W}^{-1}&0\\\
0&W\end{bmatrix}-\begin{bmatrix}A_{W}^{-1}B^{T}\\\
-I\end{bmatrix}S_{W}^{-1}\begin{bmatrix}BA_{W}^{-1}&-I\end{bmatrix}.$

Because the subtracted term is positive semidefinite, we conclude that the
eigenvalues of $\mathcal{K}^{-1}$ are less than or equal to the eigenvalues of

$\begin{bmatrix}A_{W}^{-1}&0\\\ 0&W\end{bmatrix}.$

The stated result follows. ∎

##  3 Augmentation-based bounds when $W=\gamma I$

As in section 2, we consider the augmented matrix $\mathcal{K}(W)$, but in
this case we restrict ourselves to the case where

$W=\gamma I,$

as is done in [4, 6]. For simplicity we write

$A_{\gamma}=A+\gamma B^{T}B;\qquad\mathcal{K}_{\gamma}=\mathcal{K}(\gamma I).$

In this case, the lower bound on positive eigenvalues presented in Theorem 3
reduces to $\min\left\\{\mu_{\min}(A_{\gamma}),\frac{1}{\gamma}\right\\}$.

We first consider the special case where $\operatorname{rank}(A)=n-m$ and
$\mathcal{K}$ is nonsingular. We say here that $A$ is lowest-rank because if
its rank were any lower then $\mathcal{K}$ would necessarily be singular. It
was shown in [3, 4] that $A_{\gamma}$ and $\mathcal{K}_{\gamma}$ have unique
properties, which we will use here to refine the bound on lower positive
eigenvalues given in Theorem 3. We return in Section 3.2 to the general case,
where $A$ is assumed to be rank-deficient but not lowest rank.

###  3.1 Bounds when $\operatorname{rank}(A)=n-m$

######  Theorem 4.

When $\operatorname{rank}(A)=n-m$, we have

$\mu_{\min}(A_{\gamma})\geq\rho\cdot\min\left\\{\mu_{\min}^{+}(A),\gamma\sigma_{\min}^{2}(B)\right\\},$
(7)

where $\rho\leq 1$ is a constant that does not depend on $\gamma$.

###### Proof.

We begin by writing a decomposition of $A_{\gamma}$ as was done in [4]. Let

$A=U\Lambda U^{T},\ \ B=QSV^{T}$

be the reduced (economy-size) singular value decompositions of $A$ and $B$.

The matrices $\Lambda\in\mathbb{R}^{(n-m)\times(n-m)}$ and
$U\in\mathbb{R}^{n\times(n-m)}$ comprise the eigenpairs of $A$ that correspond
to its nonzero eigenvalues, and the columns of $V\in\mathbb{R}^{n\times m}$
are the set of eigenvectors of $B^{T}B$ that correspond to its nonzero
eigenvalues. We can then write

$A_{\gamma}=P\Sigma P^{T},$ (8)

where

$P=\begin{bmatrix}U&V\end{bmatrix},\ \ \Sigma=\begin{bmatrix}\Lambda&0\\\
0&\gamma S^{2}\end{bmatrix}.$

The decomposition in (8) resembles an eigenvalue decomposition, but is not an
eigenvalue decomposition in general because the columns of $V$ will not be
orthogonal to those of $U$.

We then derive a lower bound on the eigenvalues of $A_{\gamma}$ by obtaining
an upper bound on the eigenvalues of $A_{\gamma}^{-1}$. We can write

$\mu_{\max}(A_{\gamma}^{-1})=||A_{\gamma}^{-1}||=||P^{-T}\Sigma^{-1}P^{-1}||\leq||\Sigma^{-1}||\cdot||P^{-1}||^{2}.$
(9)

The largest eigenvalue of $\Sigma^{-1}$ is equal to
$\max\left\\{\frac{1}{\mu_{\min}^{+}},\frac{1}{\gamma\sigma_{\min}^{2}}\right\\}$.
The stated result follows by setting $\rho=||P^{-1}||^{-2}$ in (9). We claim
that $\rho\leq 1$, with equality when $U$ and $V$ are mutually orthogonal
(that is, when the range of $A$ is orthogonal to the range of $B^{T}$). To
show this is the case, consider $x\in\ker(A)$. We then have

$P^{T}x=\begin{bmatrix}U^{T}x\\\ V^{T}x\end{bmatrix}=\begin{bmatrix}0\\\
V^{T}x\end{bmatrix}.$

Defining $q=P^{T}x$, since $V$ is orthogonal we have

$\|q\|\leq\|P^{-T}q\|\leq\|P^{-T}\|\|q\|,$

meaning that $||P^{-T}||$ (and therefore $||P^{-1}||$) is greater than or
equal to 1. Thus, $\rho\leq 1$. ∎

We now provide a value for $\rho=||P^{-1}||^{-2}$ in terms of the principal
angles between $\operatorname{range}(A)$ and $\operatorname{range}(B^{T})$.
Let

$\theta_{i},i=1,\ldots,m$

denote these angles. The cosines $\cos(\theta_{i})$ of these angles are given
by the singular values of $U^{T}V$ (or $V^{T}U$).

######  Lemma 5.

Let $\theta_{\min}$ denote the minimum principal angle between
$\operatorname{range}(A)$ and $\operatorname{range}(B^{T})$. Then

$||P^{-1}||=\frac{1}{\sqrt{1-\cos{(\theta_{\min})}}},$

which implies that $\rho$ defined in (7) is given by

$\rho=1-\cos{(\theta_{\min})}.$

###### Proof.

We proceed by analyzing the eigenvalues of $P^{T}P$, using the fact that

$||P^{-1}||=\frac{1}{\sqrt{\mu_{\min}({P^{T}P})}}.$

We write $P^{T}P$ in block form:

$P^{T}P=\begin{bmatrix}U^{T}\\\
V^{T}\end{bmatrix}\begin{bmatrix}U&V\end{bmatrix}=\begin{bmatrix}I&U^{T}V\\\
V^{T}U&I\end{bmatrix}.$

The (1,1)-block of $P^{T}P$ is size $(n-m)\times(n-m)$ and the (2,2)-block is
size $m\times m$. We now assume without loss of generality that $n-m\geq m$.
(If $n-m<m$, we can reorder the blocks of $P^{T}P$ such that the (1,1)-block
is larger, and use the same analysis as below.)

Letting $v=\begin{bmatrix}x^{T}&y^{T}\end{bmatrix}^{T}$ be an appropriately
partitioned eigenvector, we write the eigenvalue equations for $P^{T}P$:

$\displaystyle x+U^{T}Vy$ $\displaystyle=\lambda x;$ (10a) $\displaystyle
V^{T}Ux+y$ $\displaystyle=\lambda y.$ (10b)

There is an eigenvalue $\lambda=1$ with multiplicity $n-2m$, which we observe
by choosing $x\in\ker(V^{T}U)$ and $y=0$. For the remaining $2m$ eigenvalues,
we assume $\lambda\neq 1$. From (10a) we have $x=\frac{1}{\lambda-1}U^{T}Vy$,
which we substitute into (10b) to obtain

$y=\frac{1}{(\lambda-1)^{2}}V^{T}UU^{T}Vy.$ (11)

The eigenvalues of $V^{T}UU^{T}V$ are given by $\cos^{2}(\theta_{i})$, where
$\theta_{i}$ are the principal angles between $\operatorname{range}(A)$ and
$\operatorname{range}(B^{T})$. Thus, for each $\theta_{i}$ we can write (11)
as

$y=\frac{\cos^{2}(\theta_{i})}{(\lambda_{i}-1)^{2}}y,$

implying that

$\lambda_{i}=1\pm\cos(\theta_{i}).$

Thus each $\theta_{i}$ yields two distinct eigenvalues. Together with the
$n-2m$ eigenvalues with $\lambda=1$, this accounts for all $n$ eigenvalues of
$P^{T}P$. Therefore, the smallest eigenvalue of $P^{T}P$ is given by
$1-\cos(\theta_{\min})$; the stated result follows. ∎

We can use the results we have established for matrices with lowest-rank $A$
to derive a lower bound on the positive eigenvalues of $\mathcal{K}$ that does
not require us to know the eigenvalues of $A_{\gamma}$. We saw in Theorem 3
that for $W=\gamma I$, the bound is given by
$\min\left\\{\mu_{\min}(A_{\gamma}),\frac{1}{\gamma}\right\\}$. As $\gamma$
decreases, the value of $\mu_{\min}(A_{\gamma})$ approaches zero (because
$A_{\gamma}$ approaches $A$); thus, we achieve the best possible lower bound
when

$\frac{1}{\gamma}=\mu_{\min}(A_{\gamma}).$

Since we do not generally know the value of $\mu_{\min}(A_{\gamma})$, we can
instead select $\frac{1}{\gamma}$ to be equal to the reciprocal of the lower
bound on $\mu_{\min}(A_{\gamma})$ given by Theorem 4 and Lemma 5. That is, we
find a $\gamma$ that satisfies

$\frac{1}{\gamma}=\left(1-\cos(\theta_{\min})\right)\min\left\\{\mu_{\min}^{+},\gamma\sigma_{\min}^{2}\right\\}.$

Depending on which of the arguments to the $\min$ function is smaller, we
either have

$\frac{1}{\gamma}=\mu_{\min}^{+}(1-\cos(\theta_{\min}))$

or we have
$\frac{1}{\gamma}=(1-\cos(\theta_{\min}))\cdot\gamma\sigma_{\min}^{2}$, which
implies that

$\frac{1}{\gamma}=\sigma_{\min}\sqrt{1-\cos(\theta_{\min})}.$

Therefore, if we select

$\frac{1}{\gamma}=\min\left\\{\mu_{min}^{+}(1-\cos(\theta_{\min})),\sigma_{\min}\sqrt{1-\cos(\theta_{\min})}\right\\},$

we know that $\mu_{\min}(A_{\gamma})$ will be greater than or equal to this
value of $\frac{1}{\gamma}$. This gives the following result:

######  Theorem 6.

When $\operatorname{rank}(A)=n-m$, the positive eigenvalues of $\mathcal{K}$
are greater than or equal to

$\min\left\\{\mu_{\min}^{+}(1-\cos(\theta_{\min})),\sigma_{\min}\sqrt{1-\cos(\theta_{\min})}\right\\}.$

In some cases, more may be known about the null spaces of $A$ and $B$ than the
ranges of $A$ and $B^{T}$. For these settings, it is convenient to re-frame
the result of Theorem 6 to rely on the angle between kernels rather than the
angle between ranges. Because $\ker(A)$ and $\ker(B)$ are respectively
orthogonal to $\operatorname{range}(A)$ and $\operatorname{range}(B^{T})$, the
principal angles are the same between both pairs of subspaces. The following
result then holds.

######  Corollary 7.

Let $\operatorname{rank}(A)=n-m$ and let $\psi_{\min}$ denote the minimum
principal angle between $\ker(A)$ and $\ker(B)$. The positive eigenvalues of
$\mathcal{K}$ are greater than or equal to

$\min\left\\{\mu_{\min}^{+}(1-\cos(\psi_{\min})),\sigma_{\min}\sqrt{1-\cos(\psi_{\min})}\right\\}.$

###  3.2 Bounds when $\operatorname{rank}(A)\geq n-m$

We now return to the case in which $A$ is rank-deficient but not lowest rank,
and discuss how the results of the previous section can be extended to this
case. Let us denote the eigenvalue decomposition of $A$ by:

$A=U\Lambda U^{T}.$

Let $\Lambda^{\max}_{n-m}$ be a diagonal matrix of the $n-m$ largest
eigenvalues of $\Lambda$ and $\Lambda^{\min}_{m}$ be a diagonal matrix of the
$m$ smallest. Similarly, let $U^{\max}_{n-m}$ denote the eigenvectors
corresponding to the $n-m$ largest eigenvalues and $U^{\min}_{m}$ the
eigenvectors corresponding to the $m$ smallest eigenvalues. We then have

$A=\begin{bmatrix}U^{\max}_{n-m}&U^{\min}_{m}\end{bmatrix}\begin{bmatrix}\Lambda^{\max}_{n-m}&0\\\
0&\Lambda^{\min}_{m}\end{bmatrix}\begin{bmatrix}\left(U^{\max}_{n-m}\right)^{T}\\\
\left(U^{\min}_{m}\right)^{T}\end{bmatrix}.$ (12)

As before, if we consider a weight matrix $W=\gamma I$, a lower bound on the
positive eigenvalues of $\mathcal{K}$ is given by

$\min\left\\{\frac{1}{\gamma},\mu_{\min}(A_{\gamma})\right\\},$

as this bound does not depend on the nullity of $A$. When $A$ is not lowest
rank, the bound of Theorem 4 is not immediately applicable. However, we note
from (12) that

$A=A^{\max}_{n-m}+A^{\min}_{m},$ (13)

where
$A^{\max}_{n-m}=U^{\max}_{n-m}\Lambda^{\max}_{n-m}\left(U^{\max}_{n-m}\right)^{T}$
is semidefinite matrix with rank $n-m$ and
$A^{\min}_{m}=U^{\min}_{m}\Lambda^{\min}_{m}\left(U^{\min}_{m}\right)^{T}$ is
a semidefinite matrix with rank less than or equal to $m$. Thus, the
eigenvalues of $A_{\gamma}$ are all greater than or equal to those of

$A^{\max}_{n-m}+\gamma B^{T}B=:A^{\max}_{\gamma}.$

The eigenvalue $\mu_{n-m}$ is the smallest eigenvalue in
$\Lambda^{\max}_{n-m}$ and therefore the smallest positive eigenvalue of
$A^{\max}_{n-m}$. Let $\tilde{\theta}_{\min}$ denote the minimum principal
angle between $\operatorname{range}(A^{\max}_{n-m})$ and
$\operatorname{range}({B^{T}})$. By Theorem 4 and Lemma 5, we have

$\mu_{\min}(A_{\gamma})\geq\mu_{\min}(A^{\max}_{\gamma})\geq\left(1-\cos(\tilde{\theta}_{\min})\right)\cdot\min\left\\{\mu_{n-m},\gamma\sigma_{\min}^{2}\right\\}.$

As we did before, we can select $\frac{1}{\gamma}$ to be equal to the smaller
of these two values to obtain a lower bound on the positive eigenvalues of
$\mathcal{K}$ that does not require forming an augmented matrix. The proof of
the following theorem is similar to that of Theorem 6 and is omitted.

######  Theorem 8.

Let $A$ be semidefinite with $n-m\leq\operatorname{rank}(A)\leq n$. The
positive eigenvalues of $\mathcal{K}$ are greater than or equal to

$\min\left\\{\mu_{n-m}(1-\cos(\tilde{\theta}_{\min})),\sigma_{\min}\sqrt{1-\cos(\tilde{\theta}_{\min})}\right\\},$

where $\mu_{n-m}$ denotes the $(n-m)$-th largest eigenvalue of $A$ and
$\tilde{\theta}_{\min}$ the smallest principal angle between
$\operatorname{range}(B^{T})$ and the subspace spanned by the eigenvectors
corresponding to the $n-m$ largest eigenvalues of $A$. (Or, equivalently,
$\tilde{\theta}_{\min}$ is the smallest principal angle between $\ker(B)$ and
the subspace spanned by the eigenvectors corresponding to the $m$ smallest
eigenvalues of $A$ – see Corollary 7.)

#### Remark.

Our approach in deriving the previous result was to convert a general rank-
deficient $A$ into a lowest-rank $\tilde{A}$ by removing the part of the
spectrum corresponding to the $m$ smallest eigenvalues. However, removing this
part of the spectrum of $A$ is not always a good choice, in that it may lead
to an overly pessimistic bound. For example, consider the matrix (with $n=3$
and $m=2$):

$\mathcal{K}=\left[\begin{array}[]{c c c | c c}1&0&0&0&1\\\ 0&\alpha&0&0&0\\\ 0&0&0&1&0\\\ \hline\cr 0&0&1&0&0\\\ 1&0&0&0&0\end{array}\right]=:\begin{bmatrix}A&B^{T}\\\ B&0\end{bmatrix},$

where $0<\alpha<1$. The positive eigenvalues of $\mathcal{K}$ are $\alpha,1$,
and $\frac{1+\sqrt{5}}{2}$. The “non-removed” eigenvector $U^{\max}_{n-m}$,
which is in this case the eigenvector corresponding to $\lambda=1$, is:

$U^{\max}_{n-m}=\begin{bmatrix}1\\\ 0\\\ 0\end{bmatrix}.$

Because this eigenvector is in the range of $B^{T}$, the value
$\tilde{\theta}_{\min}$ is 0, meaning that Theorem 8 gives a bound of 0. We
would obtain a better bound if, instead of keeping the part of the spectrum of
$A$ that corresponds to the eigenvalue $\lambda=1$, we kept the portion of the
spectrum corresponding to $\lambda=\alpha$ (this would in fact give a tight
bound of $\alpha$). However, the issue of optimizing what subspace of
$\operatorname{range}(A)$ to use in order to obtain a bound is beyond the
scope of this work.

##  4 Numerical experiments

We test our eigenvalue bounds on two problems. The first is an
electromagnetics model problem described in [7]. Consider the time-harmonic
Maxwell equations in lossless media with perfectly conducting boundaries and
constant coefficients. The problem is to find the vector field $u$ and
multiplier $p$ such that

$\displaystyle\nabla\times\nabla\times u+\nabla p$ $\displaystyle=f\textrm{ in
}\Omega,$ $\displaystyle\nabla\cdot u$ $\displaystyle=0\textrm{ in }\Omega,$
$\displaystyle u\times n$ $\displaystyle=0\textrm{ on }\partial\Omega,$
$\displaystyle p$ $\displaystyle=0\textrm{ on }\partial\Omega.$

Discretizing with Nédéléc finite elements for $u$ and nodal elements for $p$
[9] yields a linear system of the form

$\begin{bmatrix}A&B^{T}\\\ B&0\end{bmatrix}\begin{bmatrix}u\\\
p\end{bmatrix}=\begin{bmatrix}g\\\ 0\end{bmatrix},$

where $A$ is a discrete curl-curl operator, $B$ is a discrete divergence
operator, and $M$ is the finite element mass matrix.

In the above-described problem, $A$ has rank $n-m$, and hence it is lowest
rank per the terminology we use in this paper. Figure 1 shows the predicted
bound (as a solid line), the actual smallest positive eigenvalue (dashed line)
for various values of $\gamma$ for a Maxwell matrix with $n=6,080$ and
$m=1,985$.

The second problem describes linear systems arising from an interior point
method (IPM) solution to a quadratic program (QP); see [10] and the references
therein for a detailed description. At each iteration of the IPM, we solve a
linear system with a matrix of the form

$\mathcal{K}=\begin{bmatrix}H+X^{-1}Z&J^{T}\\\ J&0\end{bmatrix},$

where $H$ and $J$ are respectively the Hessian and Jacobian matrices for the
QP, and $X$ and $Z$ are diagonal matrices of the current primal and dual
iterates, some entries of which go to 0 as the iterations progress. Thus, the
leading block becomes progressively more ill-conditioned as the iterations
proceed.

In Figure 2 we show the results of our bounds on the first IPM iteration on
TOMLAB111Test matrices available at https://tomopt.com/tomlab/. Problem 17 for
which the saddle-point matrix $\mathcal{K}$ is numerically singular. This
problem has $n=293$ and $m=286$. For the particular matrix shown in the
experiment below (which arises in the 12th iteration of the IPM algorithm of
[8]), there are 115 “numerically zero” eigenvalues of the leading block (which
we define as those less than machine epsilon times the largest eigenvalue of
that block).

Figure 1: Comparison of predicted and actual smallest positive eigenvalue
bounds at various values of $\gamma$ for the Maxwell matrix (lowest rank)
Figure 2: Comparison of predicted and actual smallest positive eigenvalue
bounds at various values of $\gamma$ for the IPM matrix for TOMLAB QP 17

In both cases the actual smallest positive eigenvalue
$\mu_{\min}(\mathcal{K})$ occurs precisely where
$\frac{1}{\gamma}=\mu_{\min}(A_{\gamma})$. The bounds for the Maxwell matrix
are rather tight, in the sense that they are of the same order of magnitude as
the eigenvalue (we also see this with Maxwell matrices of other sizes): the
predicted eigenvalue bound is 0.0453 while the actual smallest positive
eigenvalue is 0.0611.

The bound for the TOMLAB problem is looser: the predicted bound is
$4.716\times 10^{-7}$ while the actual smallest positive eigenvalue is
$1.817\times 10^{-4}$. Recall that our approach for deriving the bound for a
matrix with $A$ that does not have the lowest rank consisted of two steps: (1)
implicitly convert the matrix to one with a lowest-rank leading block by
“dropping” part of the spectrum of $A$ corresponding to the smallest positive
eigenvalues; and (2) estimate the lower bound for the matrix with the lowest-
rank leading block using the results of Section 3.1, using the fact that this
will also be a lower bound for the original matrix. Because our bound in the
non-lowest-rank case relies on “dropping” part of the spectrum of $A$, as
discussed in Section 3.2, we might in general expect that to lead to some
looseness in the bound.

However, the dropping is not the cause of the looseness in this case of the
TOMLAB problem, as the saddle-point matrix we obtain by simply replacing $A$
with its dropped portion $A_{n-m}^{\max}$ (defined in (13)) has almost the
same smallest positive eigenvalue as the original matrix ($1.810\times
10^{-4}$, compared with $1.817\times 10^{-4}$). Thus, the looseness in this
bound does not come from the dropping part of the spectrum of $A$ to create a
lowest-rank matrix, but rather in the estimation of the lower positive
eigenvalue bound of the modified matrix.

##  5 Conclusions

We have described a novel framework for bounding eigenvalues of saddle-point
matrices by strategically augmenting some of their blocks. We used this
approach to derive (nonzero) bounds on the lower positive eigenvalues of
saddle-point matrices with singular leading blocks. By making certain
assumptions on the augmentation parameters, we were able to derive an
eigenvalue bound that does not require the formation of an augmented matrix.

Future work may include improving the bound in the non-lowest-rank case (for
instance, by judiciously selecting the portion of the spectrum of $A$ that is
“dropped”) and using this framework to analyze the convergence of
preconditioned iterative solvers.

## References

  * [1] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numerica, 14:1–137, 2005.
  * [2] S. Bradley and C. Greif.  Eigenvalue bounds for double saddle-point systems.  https://arxiv.org/abs/2110.13328, 2021.
  * [3] R. Estrin and C. Greif.  On nonsingular saddle-point systems with a maximally rank deficient leading block.  SIAM Journal on Matrix Analysis and Applications, 36(2):367–384, 2015.
  * [4] R. Estrin and C. Greif.  Towards an optimal condition number of certain augmented Lagrangian-type saddle-point matrices.  Numerical Linear Algebra with Applications, 23(4):693–705, 2016\.
  * [5] R. Fletcher.  An Ideal Penalty Function for Constrained Optimization.  IMA Journal of Applied Mathematics, 15(3):319–342, 06 1975.
  * [6] G. H. Golub and C. Greif.  On solving block-structured indefinite linear systems.  SIAM J. Sci. Comput., 24(6):2076–2092, 2003.
  * [7] C. Greif and D. Schötzau.  Preconditioners for the discretized time-harmonic Maxwell equations in mixed form.  Numer. Linear Algebra Appl., 14(4):281–297, 2007.
  * [8] S. Mehrotra.  On the implementation of a primal-dual interior point method.  SIAM Journal on Optimization, 2(4):575–601, 1992.
  * [9] P. Monk.  Finite element methods for Maxwell’s equations.  Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2003.
  * [10] J. Nocedal and S. J. Wright.  Numerical optimization.  Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006.
  * [11] D. Ruiz, A. Sartenaer, and C. Tannier.  Refining the lower bound on the positive eigenvalues of saddle point matrices with insights on the interactions between the blocks.  SIAM Journal on Matrix Analysis and Applications, 39(2):712–736, 2018.
  * [12] T. Rusten and R. Winther.  A preconditioned iterative method for saddlepoint problems.  SIAM J. Matrix Anal. Appl., 13:887–904, 1992.