freealg.AlgebraicForm#

class freealg.AlgebraicForm(A, support=None, n=None, ratio=None, log=False, dtype='complex128', stieltjes_opt={}, inv_stieltjes_opt={}, supp_opt={})#

Algebraic surrogate for ensemble models.

Parameters:
Anumpy.ndarray

The 2D symmetric \(\mathbf{A}\). The eigenvalues of this will be computed upon calling this class. If a 1D array provided, it is assumed to be the eigenvalues of \(\mathbf{A}\).

supporttuple, default=None

The support of the density of \(\mathbf{A}\). If None, it is estimated from the minimum and maximum of the eigenvalues.

nint, default=None

When A is the eigenvalues of a matrix (not the matrix itself), and if the number of these eigenvalues are more than the intended matrix size (such as then the eigenvalues are over-sampled from multiple submatrix sizes), then the user should provide the actual size of the intended matrix using this argument.

ratiofloat, default=None

If A is a Gram (covariance) matrix that is formed by \(\mathbf{A} = \mathbf{X} \mathbf{X}^{\intercal}\) with the rectangular matrix \(\mathbb{X} \in \mathbb{R}^{p \times n}\), then this argument is the aspect ratio \(c = p / n \in (0, \infty)\). The argument is only used for deformed decompression.

logbool, default=False

If True, it is assumed the spectral density is positive-definite and its support is best represented in logarithmic scale.

dtype{'complex128', 'complex256'}, default = 'complex128'

Data type for inner computations of complex variables:

  • 'complex128': 128-bit complex numbers, equivalent of two double precision floating point.

  • 'complex256': 256-bit complex numbers, equivalent of two long double precision floating point. This option is only available on Linux machines.

stieltjes_optdict, default={}

Dictionary of options to configure the computation of Stieltjes transform from the fitted polynomial. These options are passed to StieltjesPoly class.

inv_stieltjes_optdict, default={}

Dictionary of options for computing density using Plemelj formula ( inverse Stieltjes). These options are passed to inverse_stieltjes function. The dictionary include (but not limited to) these keys:

  • deltafloat, default=1e-5

    Imaginary offset into the upper half-plane used to evaluate the Stieltjes transform for density recovery by the inverse Stieltjes transform (Plemelj formula).

  • delta_ladder_ratiofloat, default=2.0

    Geometric ratio of the imaginary offsets used for multi-level density recovery. The offsets are defined by \(\delta_j = \delta r^j\), where \(r\) is this argument.

  • delta_ladder_sizeint, default=1

    Number of imaginary offsets in the geometric ladder used for multi-level density recovery.

supp_optdict, default={}

Dictionary of parameters to pass to supp() function when support=None.

Examples

Generate a matrix:

Before working with the freealg.AlgebraicForm, we generate a matrix. This matrix can be from your empirical data. Here, we generate it from Levy distribution (see freealg.distributions.FreeLevy):

>>> # Create a distribution
>>> from freealg.distributions import FreeLevy
>>> fl = FreeLevy(t=[2.0, 5.5], w=[0.75, 1-0.75], lam=0.1, a=0,
...               sigma=0.0)

>>> # Generate a matrix realization of the distribution
>>> A = fl.matrix(size=6000, seed=0)

>>> # Sample from this matrix
>>> As = freealg.submatrix(A, size=1000, seed=0)

Create AlgebraicForm object:*

We now use this sampled matrix As to create an algebraic form. The goal is to predict the ESD of A (the larger matrix of size 6000) from the ESD of As (the smaller matrix of size 1000).

Since the matrix is generated from the Levy distribution, we now its algebraic structure \(P(z, m)\) has degrees \(d_z = \deg_z(P) = 1\) and \(d_m = \deg_m(P) = 3\) when \(\sigma = 0\).

>>> from freealg import AlgebraicForm

>>> # Fit polynomial P with degrees 1 (in z) and 3 (in m)
>>> coeffs = af.fit(deg_m=3, deg_z=1, reg=0, verbose=True)
fit residual max  : 2.3763e-01
fit residual 99.9%: 1.7992e-02

Coefficients (real)
+0.99999637 +0.95497991 +0.12108988 +0.00000000
+0.00000000 +1.00020387 +1.24161671 +0.30272471

Coefficients (imag) norm: 0.0000e+00

Stieltjes sanity check: OK

Find properties of the algebraic structure:

Once fit, we can inquiry the coefficients of the polynomial, support, branch points, etc, of the underlying distribution

>>> # Coefficients of the fitted polynomial
>>> af.coeffs.real
array([[9.999963e-01, 9.5497991e-01, 1.21089884e-01, 4.76107138e-14],
      [0.0000000e+00, 1.00020387e+00, 1.2416167e+00, 3.02724712e-01]])

# Support of the distribution of the ESD of As
>>> est_supp = af.support()
>>> print(est_supp)
[(0.019606029040223762, 1.9518480572571149)]

>>> # Atoms of distribution of the ESD of As
>>> # Output format is the tuple [(atom_location, atom_weight)]
>>> atoms = af.atoms()
>>> print(atoms)
[(-1.572739585567635e-13, 0.3999999986265035)]

>>> # Plot branch points
>>> bp = af.branch_points(tol=1e-16, real_tol=1e-16, plot=True)
>>> print(bp)
[ 1.95184806e+00+0.j          6.99218919e-01+0.08189017j
  6.99218919e-01-0.08189017j  1.96060290e-02+0.j
  -1.12675884e-13+0.j        ]

Decompress:

We now apply AlgebraicForm.decompress() to predict the ESD of the larger matrix A of size 6000. Before this, we can check if this matrix is decompressible:

>>> status = af.is_decompressible()
>>> print(status)
True

>>> # Decompress from the size 1000 of As to 6000 of A with several
>>> # intermediate sizes
>>> x = numpy.linspace(-1, 8, 300)
>>> sizes = numpy.arange(As.shape[0], A.shape[0]+1, 500)
>>> rho, x, atoms = af.decompress(
...     sizes, x=x, method='moc', min_n_times=100, newton_opt={},
...     return_atoms=True, atom_eps=1e-2, verbose=False, plot=True)
Attributes:
eignumpy.array

Eigenvalues of the matrix

nint

Initial array size (assuming a square matrix when \(\mathbf{A}\) is 2D).

ratiofloat

The aspect ratio of Gram matrix

coeffsnumpy.ndarray

The coefficients of the fitted polynomial. This attribute is available after calling fit().

supp: tuple

The predicted (or given) support \((\lambda_{\min}, \lambda_{\max})\) of the eigenvalue density.

broad_supptuple

The tuple \((\lambda_{-}, \lambda_{+})\) indicating the largest and smallest eigenvalues (considering all bulks).

est_supplist of tuples

A list of tuples (a, b) for the spectral edges of the bulks. This attribute is only available after calling support(). These spectral edges are estimated form the fitted polynomial.

_stieltjes_polycallable function

A function that computes Stieltjes transform from the fitted polynomial using homotopy method.

_momentscallable function

A function that computes evolving moments.

Methods

fit(deg_m, deg_z[, reg, r_min, r_max, n_r, ...])

Fit an algebraic structure to the input data.

support([scan_range, n_scan, refine, ...])

Estimate the spectral edges of the density.

branch_points([tol, real_tol, plot, latex, ...])

Compute global branch points.

atoms([eta, tol, real_tol, w_tol, merge_tol])

Detect atom locations and weights of distribution

density([x, eta, ac_only, plot, latex, save])

Evaluate spectral density of the fitted spectral curve.

hilbert([x, plot, latex, save])

Compute Hilbert transform of the spectral density.

stieltjes([x, y, plot, latex, save])

Compute Stieltjes transform of the spectral density

decompress(size[, x, kind, method, ...])

Free decompression of spectral density.

debug_decompress(sizes, x[, min_n_times, ...])

Plots tracked root at for point x by comparing the result from decompression and the result from all candidate roots.

candidates(size[, kind, x, eig, delta, ...])

Candidate densities of free decompression from all possible roots

is_decompressible([ratio, n_ratios, K, L, ...])

Check if the given distribution can be decompressed.

edge(t[, kind, supp, dt_max, max_iter, tol, ...])

Evolves spectral edges.

cusp(t_grid[, kind, supp, max_iter, tol, ...])

Find cusp (merge/split) point of evolving spectral edges.

eigvalsh([size, seed])

Estimate the eigenvalues.

cond([size, seed])

Estimate the condition number.

trace([size, p, seed])

Estimate the trace of a power.

slogdet([size, seed])

Estimate the sign and logarithm of the determinant.

norm([size, order, seed])

Estimate the Schatten norm.