Ctrl+K

solvers.quasisep package

solvers.quasisep package#

This subpackage implements the core linear algebra required for the QuasisepSolver, and a few extras. An intro to these matrices can be found in Matrix Computations and Semiseparable Matrices: Linear Systems.

There exist a range of definitions for quasiseparable matrices in the literature, so to be explicit, let’s select the one that we will consider in all that follows. The most suitable definition for our purposes is nearly identical to the one used by Eidelman & Gohberg (1999), with some small modifications.

Let’s start by considering an \(N \times N\) square quasiseparable matrix \(M\) with lower quasiseparable order \(m_l\) and upper quasiseparable order \(m_u\). We represent this matrix \(M\) as:

\[\begin{split}M_{ij} = \left \{ \begin{array}{ll} d_i\quad, & \mbox{if }\, i = j \\ {p}_i^T\,\left ( \prod_{k=i-1}^{j+1} {A}_k \right )\,{q}_j\quad, & \mbox{if }\, i > j \\ {h}_i^T\,\left ( \prod_{k=i+1}^{j-1} {B}_k^T \right )\,{g}_j\quad, & \mbox{if }\, i < j \\ \end{array}\right .\end{split}\]

where

  • \(i\) and \(j\) both range from \(1\) to \(N\),

  • \(d_i\) is a scalar,

  • \({p}_i\) and \({q}_j\) are both vectors with \(m_l\) elements,

  • \({A}_k\) is an \(m_l \times m_l\) matrix,

  • \({g}_j\) and \({h}_i\) are both vectors with \(m_u\) elements, and

  • \({B}_k\) is an \(m_u \times m_u\) matrix.

Comparing this definition to the one from Eidelman & Gohberg, you may notice that we have swapped the labels of \({g}_j\) and \({h}_i\), and that we’ve added an explicit transpose to \({B}_k^T\). These changes simplify the notation and implementation for symmetric matrices where, with our definition, \({g} = {p}\), \({h} = {q}\), and \({B} = {A}\).

In our definition, the product notation is a little sloppy so, to be more explicit, this is how the products expand:

\[\prod_{k=i-1}^{j+1} {A}_k \equiv {A}_{i-1}\,{A}_{i-2}\cdots{A}_{j+2}\,{A}_{j+1}\]

and

\[\prod_{k=i+1}^{j-1} {B}_k^T \equiv {B}_{i+1}^T\,{B}_{i+2}^T\cdots{B}_{j-2}^T\,{B}_{j-1}^T\]

It is often useful see a concrete example matrix. For example, under our definition, the general \(4 \times 4\) quasiseparable matrix is:

\[\begin{split}M = \left(\begin{array}{cccc} d_1 & {h}_1^T{g}_2 & {h}_1^T{B}_2^T{g}_3 & {h}_1^T{B}_2^T{B}_3^T{g}_4 \\ {p}_2^T{q}_1 & d_2 & {h}_2^T{g}_3 & {h}_2^T{B}_3^T{g}_4 \\ {p}_3^T{A}_2{q}_1 & {p}_3^T {q}_2 & d_3 & {h}_3^T{g}_4 \\ {p}_4^T{A}_3{A}_2{q}_1 & {p}_4^T{A}_3{q}_2 & {p}_4^T{q}_3 & d_4 \\ \end{array}\right)\end{split}\]

These matrices allow linear scaling for most basic linear algebra operations (that’s why we like them!), and this subpackage implements many of these building blocks.

Square Quasiseparable Matrices#

The algorithms implemented in this subpackage are are mostly based on Eidelman & Gohberg (1999) and Foreman-Mackey et al. (2017).

QSM()

The base class for all square quasiseparable matrices

DiagQSM(d)

A diagonal quasiseparable matrix

StrictLowerTriQSM(p, q, a)

A strictly lower triangular order m quasiseparable matrix

StrictUpperTriQSM(p, q, a)

A strictly upper triangular order m quasiseparable matrix

LowerTriQSM(diag, lower)

A lower triangular quasiseparable matrix

UpperTriQSM(diag, upper)

A upper triangular quasiseparable matrix

SquareQSM(diag, lower, upper)

A general square order (m1, m2) quasiseparable matrix

SymmQSM(diag, lower)

A symmetric order m quasiseparable matrix

Rectangular Quasiseparable Matrices#

While the usual definition of quasiseparable matrices is restricted to square matrices, it is useful for our purposes to also implement some algorithms for a somewhat more general class of rectangular quasiseparable matrices. These appear in the calculations for the conditional Gaussian Process when interpolating and extrapolating. We have not (yet?) worked through some of the more general operations (like scalable matrix multiplies), but those may be possible to derive.

This generalization isn’t published anywhere as far as I know (please tell me if there is a reference that you know of!), and maybe someday I’ll come up with notation that I’m satisfied with and try to write it up.

GeneralQSM(pl, ql, pu, qu, a, idx)

A rectangular (n1,n2) quasiseparable matrix with order m
Contents