Singular value expansion

sparse_ir.compute_sve(K, eps=None, cutoff=None, n_sv=None, n_gauss=None, dtype=<class 'float'>, work_dtype=None, sve_strat=None, svd_strat=None)

Perform truncated singular value expansion of a kernel.

Perform a truncated singular value expansion (SVE) of an integral kernel K : [xmin, xmax] x [ymin, ymax] -> R:

K(x, y) == sum(s[l] * u[l](x) * v[l](y) for l in (0, 1, 2, ...)),

where s[l] are the singular values, which are ordered in non-increasing fashion, u[l](x) are the left singular functions, which form an orthonormal system on [xmin, xmax], and v[l](y) are the right singular functions, which form an orthonormal system on [ymin, ymax].

The SVE is mapped onto the singular value decomposition (SVD) of a matrix by expanding the kernel in piecewise Legendre polynomials (by default by using a collocation).

Parameters:
  • K (kernel.AbstractKernel) – Integral kernel to take SVE from

  • eps (float) – Accuracy target for the basis: attempt to have singular values down to a relative magnitude of eps, and have each singular value and singular vector be accurate to eps. A work_dtype with a machine epsilon of eps**2 or lower is required to satisfy this. Defaults to 2.2e-16 if xprec is available, and 1e-8 otherwise.

  • cutoff (float) –

    Relative cutoff for the singular values. A work_dtype with machine epsilon of cutoff is required to satisfy this. Defaults to a small multiple of the machine epsilon.

    Note that cutoff and eps serve distinct purposes. cutoff reprsents the accuracy to which the kernel is reproduced, whereas eps is the accuracy to which the singular values and vectors are guaranteed.

  • n_sv (int) – Maximum basis size. If given, only at most the n_sv most significant singular values and associated singular functions are returned.

  • n_gauss (int) – Order of Legendre polynomials. Defaults to kernel hinted value.

  • dtype (np.dtype) – Data type of the result.

  • work_dtype (np.dtype) – Working data type. Defaults to a data type with machine epsilon of at most eps**2 and at most cutoff, or otherwise most accurate data type available.

  • sve_strat (AbstractSVE) – SVE to SVD translation strategy. Defaults to SamplingSVE, optionally wrapped inside of a CentrosymmSVE if the kernel is centrosymmetric.

  • svd_strat ('fast' or 'default' or 'accurate') – SVD solver. Defaults to fast (ID/RRQR) based solution when accuracy goals are moderate, and more accurate Jacobi-based algorithm otherwise.

Returns:

An SVEResult containing the truncated singular value expansion.

Expansion module

class sparse_ir.sve.AbstractSVE

Truncated singular value expansion (SVE) of an integral kernel.

Given an integral kernel K, this provides methods for computing its truncated singular value expansion (SVE), given by:

K(x, y) == sum(s[l] * u[l](x) * v[l](y) for l in range(L)),

where L is the truncation, u[l](x) is the l-th left singular function, s[l] is the l-th singular value, and v[l](y) is the l-th right singular function. The left and right singular functions form orthonormal systems on their respective spaces.

Computing the SVE involves introducing two sets of basis functions on the x and y axis and then translating the SVE into one or more matrices, the computing the singular value decomposition of those matrices, and finally postprocessing the data.

property matrices

SVD problems underlying the SVE.

postprocess(u, s, v, dtype=None)

Constructs the SVE result from the SVD

class sparse_ir.sve.CentrosymmSVE(K, eps, *, InnerSVE=None, **inner_args)

SVE of centrosymmetric kernel in block-diagonal (even/odd) basis.

For a centrosymmetric kernel K, i.e., a kernel satisfying: K(x, y) == K(-x, -y), one can make the following ansatz for the singular functions:

u[l](x) = ured[l](x) + sign[l] * ured[l](-x)
v[l](y) = vred[l](y) + sign[l] * ured[l](-y)

where sign[l] is either +1 or -1. This means that the singular value expansion can be block-diagonalized into an even and an odd part by (anti-)symmetrizing the kernel:

Keven = K(x, y) + K(x, -y)
Kodd  = K(x, y) - K(x, -y)

The l-th basis function, restricted to the positive interval, is then the singular function of one of these kernels. If the kernel generates a Chebyshev system [1], then even and odd basis functions alternate.

[1]: A. Karlin, Total Positivity (1968).

property matrices

SVD problems underlying the SVE.

postprocess(u, s, v, dtype)

Constructs the SVE result from the SVD

class sparse_ir.sve.SVEResult(u, s, v, K, eps=None)

Result of singular value expansion

class sparse_ir.sve.SamplingSVE(K, eps, *, n_gauss=None, dtype=<class 'float'>)

SVE to SVD translation by sampling technique [1].

Maps the singular value expansion (SVE) of a kernel K onto the singular value decomposition of a matrix A. This is achieved by chosing two sets of Gauss quadrature rules: (x, wx) and (y, wy) and approximating the integrals in the SVE equations by finite sums. This implies that the singular values of the SVE are well-approximated by the singular values of the following matrix:

A[i, j] = sqrt(wx[i]) * K(x[i], y[j]) * sqrt(wy[j])

and the values of the singular functions at the Gauss sampling points can be reconstructed from the singular vectors u and v as follows:

u[l,i] ≈ sqrt(wx[i]) u[l](x[i])
v[l,j] ≈ sqrt(wy[j]) u[l](y[j])

[1] P. Hansen, Discrete Inverse Problems, Ch. 3.1

property matrices

SVD problems underlying the SVE.

postprocess(u, s, v, dtype=None)

Constructs the SVE result from the SVD

sparse_ir.sve.compute(K, eps=None, cutoff=None, n_sv=None, n_gauss=None, dtype=<class 'float'>, work_dtype=None, sve_strat=None, svd_strat=None)

Perform truncated singular value expansion of a kernel.

Perform a truncated singular value expansion (SVE) of an integral kernel K : [xmin, xmax] x [ymin, ymax] -> R:

K(x, y) == sum(s[l] * u[l](x) * v[l](y) for l in (0, 1, 2, ...)),

where s[l] are the singular values, which are ordered in non-increasing fashion, u[l](x) are the left singular functions, which form an orthonormal system on [xmin, xmax], and v[l](y) are the right singular functions, which form an orthonormal system on [ymin, ymax].

The SVE is mapped onto the singular value decomposition (SVD) of a matrix by expanding the kernel in piecewise Legendre polynomials (by default by using a collocation).

Parameters:
  • K (kernel.AbstractKernel) – Integral kernel to take SVE from

  • eps (float) – Accuracy target for the basis: attempt to have singular values down to a relative magnitude of eps, and have each singular value and singular vector be accurate to eps. A work_dtype with a machine epsilon of eps**2 or lower is required to satisfy this. Defaults to 2.2e-16 if xprec is available, and 1e-8 otherwise.

  • cutoff (float) –

    Relative cutoff for the singular values. A work_dtype with machine epsilon of cutoff is required to satisfy this. Defaults to a small multiple of the machine epsilon.

    Note that cutoff and eps serve distinct purposes. cutoff reprsents the accuracy to which the kernel is reproduced, whereas eps is the accuracy to which the singular values and vectors are guaranteed.

  • n_sv (int) – Maximum basis size. If given, only at most the n_sv most significant singular values and associated singular functions are returned.

  • n_gauss (int) – Order of Legendre polynomials. Defaults to kernel hinted value.

  • dtype (np.dtype) – Data type of the result.

  • work_dtype (np.dtype) – Working data type. Defaults to a data type with machine epsilon of at most eps**2 and at most cutoff, or otherwise most accurate data type available.

  • sve_strat (AbstractSVE) – SVE to SVD translation strategy. Defaults to SamplingSVE, optionally wrapped inside of a CentrosymmSVE if the kernel is centrosymmetric.

  • svd_strat ('fast' or 'default' or 'accurate') – SVD solver. Defaults to fast (ID/RRQR) based solution when accuracy goals are moderate, and more accurate Jacobi-based algorithm otherwise.

Returns:

An SVEResult containing the truncated singular value expansion.

Singular value decomposition

sparse_ir.svd.compute(a_matrix, n_sv_hint=None, strategy='fast')

Compute thin/truncated singular value decomposition

Computes the thin/truncated singular value decomposition of a matrix A into U, s, V:

A == (U * s) @ V.T

Depending on the strategy, only as few as n_sv_hint most significant singular values may be returned, but applications should not rely on this behvaiour. The strategy parameter can be fast (RRQR/t-SVD), default (full SVD) or accurate (Jacobi rotation SVD).