guptri_py
¶
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¶
With NumPy ≤ 1.17, it may be necessary to set:
export NPY_DISTUTILS_APPEND_FLAGS=1
to fix a linking problem. See https://github.com/numpy/numpy/issues/12799.
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
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 toTrue
- 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 byguptri()
.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¶
- GUPTRI software homepage
- B. Kågström. Generalized Non-Hermitian Eigenvalue Problems, chpt. 8.7 Singular matrix pencils. 2000. https://epubs.siam.org/doi/abs/10.1137/1.9780898719581.ch8
- J. Demmel and B. Kågström. The generalized Schur decomposition of an arbitrary pencil A - λB. 1993. https://dl.acm.org/doi/10.1145/152613.152615
guptri_py
on GitHub- See also MatrixPencils.jl for similar functionality in Julia.
- See also the Matlab toolbox Matrix Canonical Structure.