Pseudospectra

C++ Header

C Header

Python Header

C++ implementation

C++ pseudospectra example driver

C++ ChunkedPseudospectra example driver

C++ TriangularPseudospectra example driver

C++ ChunkedTriangularPseudospectra example driver

The \(\epsilon\)-pseudospectrum of a square matrix \(A\) is the set of all shifts \(z\) such that \(\hat A - z\) is singular for some \(\hat A\) such that \(\| \hat A - A \|_2 < \epsilon\). In other words, \(z\) is in the \(\epsilon\)-pseudospectrum of \(A\) if the smallest singular value of \(A - z\) is less than \(\epsilon\).

Elemental currently supports two methods for computing pseudospectra: the first is a high-performance improvement of Shiu-Hong Lui’s triangularization followed by inverse iteration approach suggested in Computation of pseudospectra by continuation (please see Trefethen’s Computation of pseudospectra for a comprehensive review). In particular, Elemental begins by computing the Schur decomposition of the given matrix, which preserves the \(\epsilon\)-pseudospectrum, up to round-off error, and then simultaneously performs many Implicitly Restarted Arnoldi (IRA) iterations with the inverse normal matrix for each shift in a manner which communicates no more data than a standard triangular solve with many right-hand sides. Converged pseudospectral estimates are deflated after convergence.

The second approach is quite similar and, instead of reducing to triangular form, reduces to Hessenberg form and performs multi-shift triangular solves in a manner similar to Greg Henry’s The shifted Hessenberg system solve computation and Purkayastha et al.’s A parallel algorithm for the Sylvester-Observer Equation. This algorithm is not yet performance-tuned in Elemental, but should be preferred in massively-parallel situations where the Schur decomposition is considered infeasible.

Spectral portrait

The following routines return the norms of the shifted inverses over an automatically-determined 2D window in the complex plane (in the matrix invNormMap) with the specified x and y resolutions. The returned integer matrix corresponds to the number of iterations required for convergence at each shift in the 2D grid.

C++ API

Matrix<int> SpectralPortrait(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> SpectralPortrait(const AbstractDistMatrix<F> &A, AbstractDistMatrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

Matrix<int> TriangularSpectralPortrait(const Matrix<F> &U, Matrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> TriangularSpectralPortrait(const AbstractDistMatrix<F> &U, AbstractDistMatrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

Matrix<int> QuasiTriangularSpectralPortrait(const Matrix<Real> &U, Matrix<Real> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Real> psCtrl = PseudospecCtrl<Real>())

DistMatrix<int> QuasiTriangularSpectralPortrait(const AbstractDistMatrix<Real> &U, AbstractDistMatrix<Real> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Real> psCtrl = PseudospecCtrl<Real>())

Matrix<int> HessenbergSpectralPortrait(const Matrix<F> &H, Matrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> HessenbergSpectralPortrait(const AbstractDistMatrix<F> &H, AbstractDistMatrix<Base<F>> &invNormMap, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

C API

ElError ElSpectralPortrait_s(ElConstMatrix_s A, ElMatrix_s invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortrait_d(ElConstMatrix_d A, ElMatrix_d invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortrait_c(ElConstMatrix_c A, ElMatrix_c invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortrait_z(ElConstMatrix_z A, ElMatrix_z invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortraitDist_s(ElConstDistMatrix_s A, ElDistMatrix_s invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortraitDist_d(ElConstDistMatrix_d A, ElDistMatrix_d invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortraitDist_c(ElConstDistMatrix_c A, ElDistMatrix_c invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortraitDist_z(ElConstDistMatrix_z A, ElDistMatrix_z invNormMap, ElInt realSize, ElInt imagSize)

ElError ElSpectralPortraitX_s(ElConstMatrix_s A, ElMatrix_s invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralPortraitX_d(ElConstMatrix_d A, ElMatrix_d invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralPortraitX_c(ElConstMatrix_c A, ElMatrix_c invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralPortraitX_z(ElConstMatrix_z A, ElMatrix_z invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralPortraitXDist_s(ElConstDistMatrix_s A, ElDistMatrix_s invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralPortraitXDist_d(ElConstDistMatrix_d A, ElDistMatrix_d invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralPortraitXDist_c(ElConstDistMatrix_c A, ElDistMatrix_c invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralPortraitXDist_z(ElConstDistMatrix_z A, ElDistMatrix_z invNormMap, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

Python API

SpectralPortrait(A, realSize=100, imagSize=100, ctrl=None)

Spectral window

The following routines return the norms of the shifted inverses over a user-specified 2D window in the complex plane (in the matrix invNormMap) with the specified x and y resolutions. The returned integer matrix corresponds to the number of iterations required for convergence at each shift in the 2D grid.

C++ API

Matrix<int> SpectralWindow(const Matrix<F> &A, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> SpectralWindow(const AbstractDistMatrix<F> &A, AbstractDistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

Matrix<int> TriangularSpectralWindow(const Matrix<F> &U, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> TriangularSpectralWindow(const AbstractDistMatrix<F> &U, AbstractDistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

Matrix<int> QuasiTriangularSpectralWindow(const Matrix<Real> &U, Matrix<Real> &invNormMap, Complex<Real> center, Real realWidth, Real imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Real> psCtrl = PseudospecCtrl<Real>())

DistMatrix<int> QuasiTriangularSpectralWindow(const AbstractDistMatrix<Real> &U, AbstractDistMatrix<Real> &invNormMap, Complex<Real> center, Real realWidth, Real imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Real> psCtrl = PseudospecCtrl<Real>())

Matrix<int> HessenbergSpectralWindow(const Matrix<F> &H, Matrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int> HessenbergSpectralWindow(const AbstractDistMatrix<F> &H, AbstractDistMatrix<Base<F>> &invNormMap, Complex<Base<F>> center, Base<F> realWidth, Base<F> imagWidth, Int realSize, Int imagSize, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

C API

ElError ElSpectralWindow_s(ElConstMatrix_s A, ElMatrix_s invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindow_d(ElConstMatrix_d A, ElMatrix_d invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindow_c(ElConstMatrix_c A, ElMatrix_c invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindow_z(ElConstMatrix_z A, ElMatrix_z invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindowDist_s(ElConstDistMatrix_s A, ElDistMatrix_s invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindowDist_d(ElConstDistMatrix_d A, ElDistMatrix_d invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindowDist_c(ElConstDistMatrix_c A, ElDistMatrix_c invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindowDist_z(ElConstDistMatrix_z A, ElDistMatrix_z invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize)

ElError ElSpectralWindowX_s(ElConstMatrix_s A, ElMatrix_s invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralWindowX_d(ElConstMatrix_d A, ElMatrix_d invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralWindowX_c(ElConstMatrix_c A, ElMatrix_c invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralWindowX_z(ElConstMatrix_z A, ElMatrix_z invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralWindowXDist_s(ElConstDistMatrix_s A, ElDistMatrix_s invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralWindowXDist_d(ElConstDistMatrix_d A, ElDistMatrix_d invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

ElError ElSpectralWindowXDist_c(ElConstDistMatrix_c A, ElDistMatrix_c invNormMap, complex_float center, float realWidth, float imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_s ctrl)

ElError ElSpectralWindowXDist_z(ElConstDistMatrix_z A, ElDistMatrix_z invNormMap, complex_double center, double realWidth, double imagWidth, ElInt realSize, ElInt imagSize, ElPseudospecCtrl_d ctrl)

Python API

SpectralWindow(A, center=0, realWidth=1, imagWidth=1, realSize=100, imagSize=100, ctrl=None)

Spectral cloud

The following routines return the norms of the shifted inverses in the vector invNorms for a given set of shifts. The returned integer vector is a list of the number of iterations required for convergence of each shift.

C++ API

Matrix<int> SpectralCloud(const Matrix<F> &A, const Matrix<Complex<Base<F>>> &shifts, Matrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int, VR, STAR> SpectralCloud(const AbstractDistMatrix<F> &A, const AbstractDistMatrix<Complex<Base<F>>> &shifts, AbstractDistMatrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

Matrix<int> TriangularSpectralCloud(const Matrix<F> &U, const Matrix<Complex<Base<F>>> &shifts, Matrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int, VR, STAR> TriangularSpectralCloud(const AbstractDistMatrix<F> &U, const AbstractDistMatrix<Complex<Base<F>>> &shifts, AbstractDistMatrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int, VR, STAR> QuasiTriangularSpectralCloud(const AbstractDistMatrix<Real> &U, const AbstractDistMatrix<Complex<Real>> &shifts, AbstractDistMatrix<Real> &invNorms, PseudospecCtrl<Real> psCtrl = PseudospecCtrl<Real>())

Matrix<int> HessenbergSpectralCloud(const Matrix<F> &H, const Matrix<Complex<Base<F>>> &shifts, Matrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

DistMatrix<int, VR, STAR> HessenbergSpectralCloud(const AbstractDistMatrix<F> &H, const AbstractDistMatrix<Complex<Base<F>>> &shifts, AbstractDistMatrix<Base<F>> &invNorms, PseudospecCtrl<Base<F>> psCtrl = PseudospecCtrl<Base<F>>())

C API

ElError ElSpectralCloud_s(ElConstMatrix_s A, ElConstMatrix_c shifts, ElMatrix_s invNormMap)

ElError ElSpectralCloud_d(ElConstMatrix_d A, ElConstMatrix_z shifts, ElMatrix_d invNormMap)

ElError ElSpectralCloud_c(ElConstMatrix_c A, ElConstMatrix_c shifts, ElMatrix_s invNormMap)

ElError ElSpectralCloud_z(ElConstMatrix_z A, ElConstMatrix_z shifts, ElMatrix_d invNormMap)

ElError ElSpectralCloudDist_s(ElConstDistMatrix_s A, ElConstDistMatrix_c shifts, ElDistMatrix_s invNormMap)

ElError ElSpectralCloudDist_d(ElConstDistMatrix_d A, ElConstDistMatrix_z shifts, ElDistMatrix_d invNormMap)

ElError ElSpectralCloudDist_c(ElConstDistMatrix_c A, ElConstDistMatrix_c shifts, ElDistMatrix_s invNormMap)

ElError ElSpectralCloudDist_z(ElConstDistMatrix_z A, ElConstDistMatrix_z shifts, ElDistMatrix_d invNormMap)

ElError ElSpectralCloudX_s(ElConstMatrix_s A, ElConstMatrix_c shifts, ElMatrix_s invNormMap, ElPseudospecCtrl_s ctrl)

ElError ElSpectralCloudX_d(ElConstMatrix_d A, ElConstMatrix_z shifts, ElMatrix_d invNormMap, ElPseudospecCtrl_d ctrl)

ElError ElSpectralCloudX_c(ElConstMatrix_c A, ElConstMatrix_c shifts, ElMatrix_s invNormMap, ElPseudospecCtrl_s ctrl)

ElError ElSpectralCloudX_z(ElConstMatrix_z A, ElConstMatrix_z shifts, ElMatrix_d invNormMap, ElPseudospecCtrl_d ctrl)

ElError ElSpectralCloudXDist_s(ElConstDistMatrix_s A, ElConstDistMatrix_c shifts, ElDistMatrix_s invNormMap, ElPseudospecCtrl_s ctrl)

ElError ElSpectralCloudXDist_d(ElConstDistMatrix_d A, ElConstDistMatrix_z shifts, ElDistMatrix_d invNormMap, ElPseudospecCtrl_d ctrl)

ElError ElSpectralCloudXDist_c(ElConstDistMatrix_c A, ElConstDistMatrix_c shifts, ElDistMatrix_s invNormMap, ElPseudospecCtrl_s ctrl)

ElError ElSpectralCloudXDist_z(ElConstDistMatrix_z A, ElConstDistMatrix_z shifts, ElDistMatrix_d invNormMap, ElPseudospecCtrl_d ctrl)

Python API

SpectralCloud(A, shifts, ctrl=None)

Control structures

C++ API

class SnapshotCtrl

Int realSize

Int imagSize

Int imgSaveFreq

Int numSaveFreq

Int imgDispFreq

Negative if no snapshots should be saved/displayed, zero if only a final snapshot should be saved/displayed, and equal to \(n > 0\) if, in addition to a final snapshot, the partial results should be output roughly overy n iterations (there is no output in the middle of Impliclty Restarted Arnoldi cycles).

Int imgSaveCount

Int numSaveCount

Int imgDispCount

std::string imgBase

std::string numBase

FileFormat imgFormat

FileFormat numFormat

SnapshotCtrl()

All counters and dimensions are initially zero, all save/display “frequencies” are set to -1 (no output), the basename strings are initialized to “ps”, the image format to PNG, and the numerical format to ASCII_MATLAB.

void ResetCounts()

Resets all counters to zero

void Iterate()

Increments all counters by one

template<>
class PseudospecCtrl<Real>

bool forceComplexSchur

bool forceComplexPs

SchurCtrl<Real> schurCtrl

Int maxIts

Real tol

bool deflate

bool arnoldi

Int basisSize

bool reorthog

bool progress

SnapshotCtrl snapCtrl

template<>
class PseudospecCtrl<Base<F>>

A particular case where the datatype is the base of the potentially complex type F.

C API

ElSnapshotCtrl

ElInt realSize

ElInt imagSize

ElInt imgSaveFreq

ElInt numSaveFreq

ElInt imgDispFreq

Negative if no snapshots should be saved/displayed, zero if only a final snapshot should be saved/displayed, and equal to \(n > 0\) if, in addition to a final snapshot, the partial results should be output roughly overy n iterations (there is no output in the middle of Impliclty Restarted Arnoldi cycles).

ElInt imgSaveCount

ElInt numSaveCount

ElInt imgDispCount

const char* imgBase

const char* numBase

ElFileFormat imgFormat

ElFileFormat numFormat

ElError ElSnapshotCtrlDefault(ElSnapshotCtrl* ctrl)

ElError ElSnapshotCtrlDestroy(ElSnapshotCtrl* ctrl)

ElPseudospecCtrl_s

bool forceComplexSchur

bool forceComplexPs

ElSchurCtrl_s schurCtrl

ElInt maxIts

float tol

bool deflate

bool arnoldi

ElInt basisSize

bool reorthog

bool progress

ElSnapshotCtrl snapCtrl

ElPseudospecCtrl_d

bool forceComplexSchur

bool forceComplexPs

ElSchurCtrl_s schurCtrl

ElInt maxIts

double tol

bool deflate

bool arnoldi

ElInt basisSize

bool reorthog

bool progress

ElSnapshotCtrl snapCtrl

ElError ElPseudospecCtrlDefault_s(ElPseudospecCtrl_s* ctrl)

ElError ElPseudospecCtrlDefault_d(ElPseudospecCtrl_d* ctrl)

ElError ElPseudospecCtrlDestroy_s(ElPseudospecCtrl_s* ctrl)

ElError ElPseudospecCtrlDestroy_d(ElPseudospecCtrl_d* ctrl)

Python API

TODO