guptri_py

Documentation Status

A GUPTRI wrapper for NumPy and SageMath

This Python package provides Python bindings for the software GUPTRI and can be used with NumPy and, optionally, SageMath.

GUPTRI is a Fortran library by Jim Demmel and Bo Kågström for robust computation of generalized eigenvalues of singular matrix pencils. Standard tools like LAPACK do not reliably handle singular generalized eigenvalue problems.

GUPTRI solves this by computing a generalized block upper triangular form (generalized Schur staircase form) of a matrix pencil, revealing the Kronecker structure of the pencil. For details, see the documentation and the references therein.

Examples

See the examples and documentation at https://guptri-py.readthedocs.io.

Installation

Requirements: NumPy and, optionally, SageMath (tested with Sage 9.6 on Arch Linux and with earlier versions on macOS).

First, clone the repository from GitHub:

git clone https://github.com/mwageringel/guptri_py.git && cd guptri_py

To install with Python 3 and NumPy, run the following command:

pip3 install --upgrade --no-index -v .

Alternatively, for use with Sage, run this command:

sage -pip install --upgrade --no-index -v .

To install into the Python user install directory (no root access required), use:

sage -pip install --upgrade --no-index -v --user .

After successful installation, run the tests with Sage:

sage -t guptri_py

Issues



Documentation for guptri_py

For a matrix pencil A - λ B, the rank-reducing λ are called eigenvalues. If the pencil is regular, these are the roots of \(\det(A - λB)\). If the determinant is constantly zero for all λ, the pencil is singular. In this case, the eigenvalues are found by computing the GUPTRI form of the pencil.

A pencil A - λ B in GUPTRI form is block upper triangular such that

\[\begin{split}A - λ B = \begin{pmatrix} A_r - λ B_r & * & * & * & * \\ 0 & A_0 - λ B_0 & * & * & * \\ 0 & 0 & A_f - λ B_f & * & * \\ 0 & 0 & 0 & A_∞ - λ B_∞ & * \\ 0 & 0 & 0 & 0 & A_l - λ B_l \end{pmatrix}\end{split}\]

where

  • \(A_r - λ B_r\) has all right singular structure,
  • \(A_0 - λ B_0\) has all Jordan structure of the 0 eigenvalue,
  • \(A_f - λ B_f\) has all Jordan structure for non-zero finite eigenvalues,
  • \(A_∞ - λ B_∞\) has all Jordan structure of infinite eigenvalue,
  • \(A_l - λ B_l\) has all left-singular structure of the pencil.
guptri_py.guptri(A, B, *, epsu=None, gap=1000, zero=True, subdivide=True, compact=True, part=None)

Compute the GUPTRI form of two matrices A and B.

Parameters:
  • A – NumPy array or Sage matrix of size m × n
  • B – NumPy array or Sage matrix of size m × n
  • epsu (float, optional) – relative uncertainty in data, defaults to about 1e-8
  • gap (float, optional) – used for rank decisions in SVDs, defaults to 1000, should be at least 1
  • zero (boolean, optional) – truncate small singular values to zero during reduction process, defaults to True
  • subdivide (boolean, optional) – when using Sage, return block matrices, defaults to True
  • compact (boolean, optional) – strip diagonal blocks of size zero from the result when using subdivide, defaults to True
  • part (list or tuple of indices 0..4, optional) –

    return only the submatrices corresponding to specific diagonal blocks, for example

    • [0] - minimal left and right reducing subspaces
    • [0,1,2,3] - maximal left and right reducing subspaces
Return type:

5-tuple of NumPy arrays or Sage matrices (depending on input)

Returns:

S, T, P, Q, kstr where

  • S - λ T is in GUPTRI form,
  • P, Q are unitary,
  • kstr contains the Kronecker structure of the pencil (see the References),

such that

\[P^* (A - λ B) Q = S - λ T.\]

The eigenvalues are the ratios of the diagonal elements of the regular part of the pencil, i.e. the three central square diagonal blocks of S and T.

The leading columns of P and Q form orthogonal bases of the corresponding left and right reducing subspaces, respectively.

EXAMPLES:

Using NumPy:

>>> import numpy as np
>>> from guptri_py import *
>>> A = np.array([[22,34,31,31,17],
...               [45,45,42,19,29],
...               [39,47,49,26,34],
...               [27,31,26,21,15],
...               [38,44,44,24,30]], float)
>>> B = np.array([[13,26,25,17,24],
...               [31,46,40,26,37],
...               [26,40,19,25,25],
...               [16,25,27,14,23],
...               [24,35,18,21,22]], float)
>>> S, T, P, Q, kstr = guptri(A, B); S, T, kstr
(array([[   0.        ,    0.        ,  -31.21794153,   69.71856621, -142.67272727],
        [   0.        ,    0.        ,    0.        ,   15.85324665,  -42.10391654],
        [   0.        ,    0.        ,    0.        ,  -15.56393985,   -3.47120411],
        [   0.        ,    0.        ,    0.        ,    0.        ,   13.59793536],
        [   0.        ,    0.        ,    0.        ,    0.        ,    0.        ]]),
 array([[   0.        ,  -21.59418044,  -22.16602894,   44.51530913, -110.98181818],
        [   0.        ,    0.        ,  -25.05806277,   19.3190495 ,  -43.89049025],
        [   0.        ,    0.        ,    0.        ,   -7.78196993,    9.20409929],
        [   0.        ,    0.        ,    0.        ,    0.        ,    0.        ],
        [   0.        ,    0.        ,    0.        ,    0.        ,    0.        ]]),
 array([[ 2,  1,  0, -1,  2,  0, -1,  1, -1],
        [ 1,  1,  0, -1,  1,  0, -1,  1, -1]], dtype=int32))
>>> np.linalg.norm(A - P.dot(S.dot(Q.T.conj()))) < 1e-12
True
>>> np.linalg.norm(B - P.dot(T.dot(Q.T.conj()))) < 1e-12
True

We extract the block sizes and conclude that there are two 0 eigenvalues, one finite eigenvalue at 2 and one infinite eigenvalue:

>>> kcf_blocks(kstr)
array([[0, 2, 1, 1, 1],
       [1, 2, 1, 1, 0]])
>>> S[2,3] / T[2,3]  # tol 1e-13
2.0

Using Sage:

sage: from guptri_py import *
sage: A = matrix(A); B = matrix(B)
sage: S, T, P, Q, kstr = guptri(A, B)
sage: S  # tol 1e-13
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0 -31.217941525933004|  69.71856620722244|-142.67272727272717|]
[                0.0|                0.0                 0.0|  15.85324664502821| -42.10391654473494|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0|-15.563939851607742|-3.4712041141556753|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0|                0.0| 13.597935363644218|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0|                0.0|                0.0|]
sage: T  # tol 1e-13
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|-21.594180436280595 -22.166028942864237|  44.51530913141558| -110.9818181818181|]
[                0.0|                0.0 -25.058062773982755|  19.31904950452907| -43.89049025352817|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0| -7.781969925803871|  9.204099294191801|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0|                0.0|                0.0|]
[-------------------+---------------------------------------+-------------------+-------------------+]
[                0.0|                0.0                 0.0|                0.0|                0.0|]
sage: P  # tol 1e-13
[|  0.2696799449852967  0.45203748616024275|  0.8072453264583416|  0.0225657325979767|-0.26604625350141314]
[|  0.4045199174779453   0.5100693255997332| -0.4984905633467818|  0.5399884462530116|-0.19003303821529563]
[|  0.6741998624632421    -0.38178841736507|  0.1704780900615978|  0.2135798916315633|  0.5700991146458833]
[| 0.13483997249264842   0.5619925503613826|-0.17710779356399256| -0.6188095967360872|  0.5016872208883796]
[|  0.5393598899705936 -0.27183335316392954| -0.1985754049050833| -0.5285466663441797|  -0.562497793117271]
sage: Q  # tol 1e-13
[ -0.4043680421515943| 0.15099829724341438   0.4636322195173679|  0.7252612584539466| -0.2696799449852965|]
[  0.6984538909891174|  -0.409711992194672   0.1200055765268443| 0.19745332690892747| -0.5393598899705937|]
[  0.1470429244187616|  0.6895752857826296   -0.563291481656615| 0.14808999518169602| -0.4045199174779448|]
[ -0.4043680421515946|-0.03098662125842025  0.18594270648664857| -0.5885628013631496| -0.6741998624632417|]
[-0.40436804215159394| -0.5769413767639264  -0.6471258326055097|  0.2582081967270594|-0.13483997249264815|]
sage: kstr
[ 2  1  0 -1  2  0 -1  1 -1]
[ 1  1  0 -1  1  0 -1  1 -1]
sage: (A - P * S * Q.H).norm() < 1e-12
True
sage: (B - P * T * Q.H).norm() < 1e-12
True

A rectangular example:

sage: A = matrix([[0,1,0], [0,0,2]])
sage: B = matrix([[0,0,0], [0,0,3]])
sage: guptri(A, B)
(
[----+----+----]  [---+---+---]
[ 0.0| 2.0| 0.0]  [0.0|3.0|0.0]              [-1.0| 0.0|-0.0]
[----+----+----]  [---+---+---]  [|0.0|1.0]  [-0.0| 0.0|-1.0]
[ 0.0| 0.0|-1.0], [0.0|0.0|0.0], [|1.0|0.0], [ 0.0| 1.0|-0.0],
<BLANKLINE>
[ 1 -1  1  0 -1  1 -1]
[ 0 -1  1  0 -1  1 -1]
)

Preserving empty diagonal blocks:

sage: S, T = guptri(A, B, compact=False)[:2]; S, T
(
[----++----+----+]  [---++---+---+]
[----++----+----+]  [---++---+---+]
[ 0.0|| 2.0| 0.0|]  [0.0||3.0|0.0|]
[----++----+----+]  [---++---+---+]
[ 0.0|| 0.0|-1.0|]  [0.0||0.0|0.0|]
[----++----+----+], [---++---+---+]
)
sage: S.subdivision(3, 3), T.subdivision(3, 3)
([-1.0], [0.0])

An example showing that SciPy, using LAPACK, is not suited for solving singular eigenvalue problems:

sage: A = matrix(RDF, [[1, 2e-16], [3e-16, 0]])
sage: B = matrix(RDF, [[1, 1e-16], [1e-16, 0]])
sage: import scipy.linalg
sage: scipy.linalg.eigvals(A, B)
array([2.+0.j, 3.+0.j])
sage: scipy.linalg.eigvals(A, B, homogeneous_eigvals=True)
array([[2.e-16+0.j, 3.e-16+0.j],
       [1.e-16+0.j, 1.e-16+0.j]])

Using GUPTRI instead, we find a single eigenvalue at 1:

sage: guptri(A, B)[:2]
(
[---+---+]  [---+---+]
[0.0|1.0|]  [0.0|1.0|]
[---+---+]  [---+---+]
[0.0|0.0|], [0.0|0.0|]
)

The minimal and maximal left and right reducing subspaces:

sage: guptri(A, B, part=[0])[2:4]
(
    [0.0]
[], [1.0]
)
sage: guptri(A, B, part=range(4))[2:4]
(
[|  1.0]  [0.0|1.0]
[|1e-16], [1.0|0.0]
)

TESTS:

sage: import guptri_py
sage: TestSuite(guptri_py.guptri_py._tests_sage()).run(skip='_test_pickling')
guptri_py.kcf_blocks(kstr)

Compute the block sizes in the Kroncker structure kstr returned by guptri().

Parameters:kstr (NumPy array or Sage matrix) – the Kronecker structure computed by guptri()
Returns:NumPy array of size 2 × 5 where index (0,k) and (1,k) are the number of rows and columns of diagonal block k, respectively

EXAMPLES:

>>> import numpy as np
>>> from guptri_py import *
>>> A = np.array([[0,1,0], [0,0,2]], float)
>>> B = np.array([[0,0,0], [0,0,3]], float)
>>> kstr = guptri(A, B)[-1]
>>> kcf_blocks(kstr)
array([[0, 0, 1, 1, 0],
       [1, 0, 1, 1, 0]])


References