API — transforms

Numerical transforms for matrix-/operator-valued free probability.

Primary API (v0.1.0)

Cauchy transforms (matrix-valued):

• :func:cauchy_matrix_semicircle — \(G(z)\) for \(S = \sum_i A_i \otimes X_i\)
• :func:cauchy_kronecker — \(G_a(w) = E[(w - a \otimes X)^{-1}]\) for any scalar distribution (fast, \(O(n^3)\))
• :func:cauchy_kronecker_semicircle — specialization to \(X\) semicircular
• :func:cauchy_biased_matrix_semicircle — \(G(z)\) for \(S = a_0 + \sum_i A_i \otimes X_i\)
• :func:cauchy_polynomial — \(G(z)\) for a self-adjoint polynomial \(p(X_1,\dots,X_s)\) via linearization

Scalar densities:

• :func:matrix_semicircle_density — density of \(\sum_i A_i \otimes X_i\)
• :func:biased_matrix_semicircle_density — density of \(a_0 + \sum_i A_i \otimes X_i\)
• :func:polynomial_density — density of \(p(X_1,\dots,X_s)\)

Scalar (classical) helpers:

• :func:semicircle_density_scalar — Wigner semicircle density
• :func:semicircle_cauchy_scalar — scalar Cauchy transform of the semicircle law

Utilities:

• :func:subordination_kronecker — subordination function \(\omega_1(b)\) for free additive convolution
• :func:lambda_eps — regularized spectral parameter \(\Lambda_\varepsilon(z)\)

cauchy_matrix_semicircle(z, A, G0=None, tol=1e-10, maxiter=500)

Matrix-valued Cauchy transform of the matrix semicircle \(\sum_i A_i \otimes X_i\).

Solves the operator-valued semicircle equation $$ z\,G \;=\; I \;+\; \eta(G)\,G, \qquad \Im z>0, $$ by fixed-point iteration using the half-averaged map $$ G \;\mapsto\; \tfrac12\Big[\,G + (\,zI - \eta(G)\,)^{-1}\Big], $$ where \(\eta(B)=\sum_{i=1}^s A_i B A_i^\ast\).

This follows the numerical damping suggested by Helton--Rashidi Far--Speicher (IMRN 2007).

Parameters:

Name	Type	Description	Default
z	complex	Spectral parameter with \(\Im z>0\).	required
A	sequence of $(n,n)$ arrays or stacked $(s,n,n)$ array	Kraus operators \(A_i\) defining \(\eta(B)=\sum_i A_i B A_i^\ast\).	required
G0	(n,n) array	Initial iterate (defaults to \(-iI\)).	None
tol	float	Relative fixed-point tolerance.	1e-10
maxiter	int	Maximum iterations.	500

Returns:

Type	Description
(n, n) ndarray	Approximate solution \(G(z)\).

Notes

The residual \(R=zG-I-\eta(G)G\) should be small at convergence.

References

• J. W. Helton, R. Rashidi Far, R. Speicher, Operator-valued Semicircular Elements, IMRN (2007).

Source code in src/free_matrix_laws/transforms.py

def cauchy_matrix_semicircle(z: complex, A, G0: np.ndarray | None = None,
                             tol: float = 1e-10, maxiter: int = 500) -> np.ndarray:
    r'''
    Matrix-valued Cauchy transform of the matrix semicircle $\sum_i A_i \otimes X_i$.

    Solves the operator-valued semicircle equation
    $$ z\,G \;=\; I \;+\; \eta(G)\,G, \qquad \Im z>0, $$
    by fixed-point iteration using the half-averaged map
    $$ G \;\mapsto\; \tfrac12\Big[\,G + (\,zI - \eta(G)\,)^{-1}\Big], $$
    where $\eta(B)=\sum_{i=1}^s A_i B A_i^\ast$.

    This follows the numerical damping suggested by Helton--Rashidi Far--Speicher (IMRN 2007).

    Parameters
    ----------
    z : complex
        Spectral parameter with $\Im z>0$.
    A : sequence of $(n,n)$ arrays or stacked $(s,n,n)$ array
        Kraus operators $A_i$ defining $\eta(B)=\sum_i A_i B A_i^\ast$.
    G0 : (n,n) array, optional
        Initial iterate (defaults to $-iI$).
    tol : float
        Relative fixed-point tolerance.
    maxiter : int
        Maximum iterations.

    Returns
    -------
    (n, n) ndarray
        Approximate solution $G(z)$.

    Notes
    -----
    The residual $R=zG-I-\eta(G)G$ should be small at convergence.

    References
    ----------
    * J. W. Helton, R. Rashidi Far, R. Speicher, *Operator-valued Semicircular
      Elements*, IMRN (2007).
    '''
    n = (A[0].shape[0] if isinstance(A, (list, tuple)) else A.shape[-1])
    G = (-1j * np.eye(n)) if G0 is None else np.array(G0, dtype=complex)

    for _ in range(maxiter):
        G_next = _hfs_map(G, z, A)
        if la.norm(G_next - G) <= tol * (1 + la.norm(G)):
            return G_next
        G = G_next
    return G

cauchy_biased_matrix_semicircle(z, a0, A, G0=None, tol=1e-12, maxiter=5000, relax=0.5, return_info=False)

Matrix-valued Cauchy transform of the biased matrix semicircle \(S = a_0 + \sum_i A_i \otimes X_i\).

Solves $$ G \;=\; (z I - a_0 - \eta(G))^{-1}, \qquad \Im z>0 , $$ via the relaxed iteration $$ G_{k+1} \;=\; (1-\text{relax})\,G_k \;+\; \text{relax}\,[\,z I - b\,\eta(G_k)\,]^{-1} b, \quad b=z(z I - a_0)^{-1}. $$

Convergence is checked via the residual $$ R(G):= (z I - a_0)\,G - I - \eta(G)\,G, $$ stopping when \(\|R(G)\|_{\mathrm{F}} \le \text{tol}\).

Parameters:

Name	Type	Description	Default
z	complex	Spectral parameter with \(\Im z>0\).	required
a0	(n,n) ndarray	Bias matrix.	required
A	sequence[(n,n)] or (s,n,n) ndarray	Kraus operators for \(\eta\).	required
G0	(n,n) ndarray	Warm start; if None, uses \((\Im z)^{-1} i\,I\).	None
tol	float	Frobenius-norm tolerance on the residual.	1e-12
maxiter	int	Iteration cap.	5000
relax	float	Averaging parameter in \((0,1]\).	0.5
return_info	bool	If True, also return a dict with residual and iterations.	False

Returns:

Name	Type	Description
G	(n,n) ndarray
info	dict (only if return_info=True)	Keys: residual, iters.

Source code in src/free_matrix_laws/transforms.py

def cauchy_biased_matrix_semicircle(
    z: complex,
    a0: np.ndarray,
    A,
    G0: np.ndarray | None = None,
    tol: float = 1e-12,
    maxiter: int = 5000,
    relax: float = 0.5,
    return_info: bool = False,
):
    r'''
    Matrix-valued Cauchy transform of the biased matrix semicircle
    $S = a_0 + \sum_i A_i \otimes X_i$.

    Solves
    $$
      G \;=\; (z I - a_0 - \eta(G))^{-1},
      \qquad \Im z>0 ,
    $$
    via the relaxed iteration
    $$
      G_{k+1} \;=\; (1-\text{relax})\,G_k \;+\; \text{relax}\,[\,z I - b\,\eta(G_k)\,]^{-1} b,
      \quad b=z(z I - a_0)^{-1}.
    $$

    Convergence is checked via the residual
    $$
      R(G):= (z I - a_0)\,G - I - \eta(G)\,G,
    $$
    stopping when $\|R(G)\|_{\mathrm{F}} \le \text{tol}$.

    Parameters
    ----------
    z : complex
        Spectral parameter with $\Im z>0$.
    a0 : (n,n) ndarray
        Bias matrix.
    A : sequence[(n,n)] or (s,n,n) ndarray
        Kraus operators for $\eta$.
    G0 : (n,n) ndarray, optional
        Warm start; if None, uses $(\Im z)^{-1} i\,I$.
    tol : float, default 1e-12
        Frobenius-norm tolerance on the residual.
    maxiter : int, default 5000
        Iteration cap.
    relax : float, default 0.5
        Averaging parameter in $(0,1]$.
    return_info : bool, default False
        If True, also return a dict with residual and iterations.

    Returns
    -------
    G : (n,n) ndarray
    info : dict (only if return_info=True)
        Keys: ``residual``, ``iters``.
    '''
    n = a0.shape[0]
    I = np.eye(n, dtype=complex)
    if G0 is None:
        eps_imag = float(np.imag(z))
        if eps_imag <= 0:
            eps_imag = 1e-2
        G = -1j * I / eps_imag
    else:
        G = G0.astype(complex, copy=True)

    for k in range(1, maxiter + 1):
        G = _hfsb_map(G, z, a0, A, relax=relax)
        r = (z * I - a0) @ G - I - eta(G, A) @ G
        res = la.norm(r, 'fro')
        if res <= tol:
            if return_info:
                return G, {"residual": res, "iters": k}
            return G
    if return_info:
        return G, {"residual": res, "iters": maxiter}
    return G

cauchy_polynomial(z, a0, A, *, eps_reg=1e-06, block_size=1, G0=None, tol=1e-10, maxiter=10000, return_info=False)

Cauchy transform of a self-adjoint polynomial \(p(X_1,\dots,X_s)\) via linearization.

Given a self-adjoint linearization $$ L_p = a_0 + \sum_{i=1}^s A_i \otimes X_i, $$ with semicircular \(X_i\), we form $$ b_\varepsilon(z) := z\big(\Lambda_\varepsilon(z)-a_0\big)^{-1}, \qquad \eta(B) := \sum_{i=1}^s A_i\,B\,A_i^\ast, $$ where $$ \Lambda_\varepsilon(z)=\operatorname{diag} \big(z I_k,\ i\varepsilon\, I_{n-k}\big), \qquad k=\text{block_size}. $$

The iteration uses the half-averaged update $$ G_{new} = \frac12\Big[G + \big(zI - b_\varepsilon(z)\,\eta(G)\big)^{-1} b_\varepsilon(z)\Big], $$ and stops when \(\|G_{new}-G\|_F \le \text{tol}\,\|G\|_F\).

Parameters:

Name	Type	Description	Default
z	complex	Spectral parameter with \(\Im z>0\) (typically \(z=x+i\,\varepsilon\)).	required
a0	(n,n) array	The bias / constant term of the linearization.	required
A	sequence of (n,n) arrays or stacked array (s,n,n)	Coefficients \(A_i\) defining \(\eta(B)=\sum_i A_i B A_i^\ast\).	required
eps_reg	float	Regularization parameter in \(\Lambda_\varepsilon(z)\) (lower block).	1e-6
block_size	int	Size \(k\) of the distinguished top-left block.	1
G0	(n,n) array	Initial iterate. If None, uses \(G_0 = (1/z)I\).	None
tol	float	Relative tolerance.	1e-10
maxiter	int	Maximum number of iterations.	10000
return_info	bool	If True, also return a dict with diagnostics.	False

Returns:

Name	Type	Description
G	(n,n) array	Approximation to the fixed point \(G(z,b_\varepsilon(z))\).
info	dict(optional)	Keys: 'iters', 'last_diff'.

Notes

After computing \(G(z,b_\varepsilon(z))\), the scalar Cauchy transform of \(p\) is obtained from the distinguished corner via $$ m_p(z) \approx \frac{1}{k}\,\mathrm{tr}\, \big(G(z,b_\varepsilon(z))\big)_{11}, $$ with \(k=\text{block\_size}\) and \((\cdot)_{11}\) the top-left \(k\times k\) block.

Source code in src/free_matrix_laws/transforms.py

def cauchy_polynomial(
    z: complex,
    a0: np.ndarray,
    A,
    *,
    eps_reg: float = 1e-6,
    block_size: int = 1,
    G0: np.ndarray | None = None,
    tol: float = 1e-10,
    maxiter: int = 10_000,
    return_info: bool = False,
):
    r'''
    Cauchy transform of a self-adjoint polynomial $p(X_1,\dots,X_s)$ via linearization.

    Given a self-adjoint linearization
    $$
      L_p = a_0 + \sum_{i=1}^s A_i \otimes X_i,
    $$
    with semicircular $X_i$, we form
    $$
      b_\varepsilon(z) := z\big(\Lambda_\varepsilon(z)-a_0\big)^{-1},
      \qquad
      \eta(B) := \sum_{i=1}^s A_i\,B\,A_i^\ast,
    $$
    where
    $$
    \Lambda_\varepsilon(z)=\operatorname{diag}
        \big(z I_k,\ i\varepsilon\, I_{n-k}\big),
        \qquad k=\text{block\_size}.
    $$

    The iteration uses the half-averaged update
    $$
      G_{new} = \frac12\Big[G + \big(zI - b_\varepsilon(z)\,\eta(G)\big)^{-1}
                 b_\varepsilon(z)\Big],
    $$
    and stops when $\|G_{new}-G\|_F \le \text{tol}\,\|G\|_F$.

    Parameters
    ----------
    z : complex
        Spectral parameter with $\Im z>0$ (typically $z=x+i\,\varepsilon$).
    a0 : (n,n) array
        The bias / constant term of the linearization.
    A : sequence of (n,n) arrays or stacked array (s,n,n)
        Coefficients $A_i$ defining $\eta(B)=\sum_i A_i B A_i^\ast$.
    eps_reg : float, default 1e-6
        Regularization parameter in $\Lambda_\varepsilon(z)$ (lower block).
    block_size : int, default 1
        Size $k$ of the distinguished top-left block.
    G0 : (n,n) array, optional
        Initial iterate. If None, uses $G_0 = (1/z)I$.
    tol : float, default 1e-10
        Relative tolerance.
    maxiter : int, default 10000
        Maximum number of iterations.
    return_info : bool, default False
        If True, also return a dict with diagnostics.

    Returns
    -------
    G : (n,n) array
        Approximation to the fixed point $G(z,b_\varepsilon(z))$.
    info : dict (optional)
        Keys: 'iters', 'last_diff'.

    Notes
    -----
    After computing $G(z,b_\varepsilon(z))$, the scalar Cauchy transform of $p$
    is obtained from the distinguished corner via
    $$
      m_p(z) \approx \frac{1}{k}\,\mathrm{tr}\,
          \big(G(z,b_\varepsilon(z))\big)_{11},
    $$
    with $k=\text{block\_size}$ and $(\cdot)_{11}$ the top-left $k\times k$ block.
    '''
    a0 = np.asarray(a0)
    if a0.ndim != 2 or a0.shape[0] != a0.shape[1]:
        raise ValueError(f"a0 must be square; got {a0.shape!r}")
    n = a0.shape[0]

    if G0 is None:
        G = (1.0 / complex(z)) * np.eye(n, dtype=complex)
    else:
        G = np.asarray(G0, dtype=complex)
        if G.shape != (n, n):
            raise ValueError(f"G0 must have shape {(n,n)!r}; got {G.shape!r}")

    last_diff = np.inf
    for k in range(1, maxiter + 1):
        G1 = _hfsc_map(G, z, a0, A, eps_reg=eps_reg, block_size=block_size)
        diff = la.norm(G1 - G, "fro")
        denom = max(1.0, la.norm(G, "fro"))
        last_diff = diff

        if diff <= tol * denom:
            G = G1
            break
        G = G1

    if return_info:
        return G, {"iters": k, "last_diff": float(last_diff)}
    return G

matrix_semicircle_density(x, A, eps=0.01, G0=None, tol=1e-10, maxiter=10000, a0=None)

Scalar density of the matrix semicircle \(\sum_i A_i \otimes X_i\) (optionally with bias \(a_0\)).

Unbiased case (\(a_0\) is None): compute \(G(z)\) for \(z=x+i\varepsilon\) from $$ z\,G \;=\; I \;+\; \eta(G)\,G, \qquad \eta(B)=\sum_{i=1}^s A_i B A_i^\ast, $$ then return the scalar density $$ f(x) \;=\; -\frac{1}{\pi}\,\Im!\left(\frac{1}{n}\, \mathrm{tr}\,G(x+i\varepsilon)\right). $$

Biased case (\(a_0 \ne 0\)): compute \(G_{a_0+X}(z)\) via :func:cauchy_biased_matrix_semicircle and apply the same inversion.

Parameters:

Name	Type	Description	Default
x	float	Real evaluation point.	required
A	sequence of $(n,n)$ arrays or stacked array $(s,n,n)$	Kraus operators \(A_i\).	required
eps	float	Imaginary offset \(\varepsilon>0\) for \(z=x+i\varepsilon\).	1e-2
G0	(n,n) array	Initial iterate (passed to the solver).	None
tol	float	Relative fixed-point tolerance.	1e-10
maxiter	int	Maximum iterations.	10000
a0	(n,n) array or ``None``	Bias matrix. If provided, computes density for \(a_0 + \sum_i A_i \otimes X_i\).	None

Returns:

Type	Description
float	Approximation to \(f(x)\).

Source code in src/free_matrix_laws/transforms.py

def matrix_semicircle_density(
    x: float,
    A,
    eps: float = 1e-2,
    G0=None,
    tol: float = 1e-10,
    maxiter: int = 10_000,
    a0=None,
) -> float:
    r'''
    Scalar density of the matrix semicircle $\sum_i A_i \otimes X_i$
    (optionally with bias $a_0$).

    Unbiased case ($a_0$ is ``None``): compute $G(z)$ for $z=x+i\varepsilon$ from
    $$ z\,G \;=\; I \;+\; \eta(G)\,G, \qquad \eta(B)=\sum_{i=1}^s A_i B A_i^\ast, $$
    then return the scalar density
    $$ f(x) \;=\; -\frac{1}{\pi}\,\Im\!\left(\frac{1}{n}\,
       \mathrm{tr}\,G(x+i\varepsilon)\right). $$

    Biased case ($a_0 \ne 0$): compute $G_{a_0+X}(z)$ via
    :func:`cauchy_biased_matrix_semicircle` and apply the same inversion.

    Parameters
    ----------
    x : float
        Real evaluation point.
    A : sequence of $(n,n)$ arrays or stacked array $(s,n,n)$
        Kraus operators $A_i$.
    eps : float, default 1e-2
        Imaginary offset $\varepsilon>0$ for $z=x+i\varepsilon$.
    G0 : (n,n) array, optional
        Initial iterate (passed to the solver).
    tol : float, default 1e-10
        Relative fixed-point tolerance.
    maxiter : int, default 10000
        Maximum iterations.
    a0 : (n,n) array or ``None``
        Bias matrix. If provided, computes density for $a_0 + \sum_i A_i \otimes X_i$.

    Returns
    -------
    float
        Approximation to $f(x)$.
    '''
    if eps <= 0:
        raise ValueError("eps must be > 0")

    if isinstance(A, np.ndarray) and A.ndim == 3:
        n = A.shape[-1]
    elif isinstance(A, np.ndarray) and A.ndim == 2:
        n = A.shape[0]
        A = A[None, ...]
    else:
        A = list(A)
        n = A[0].shape[0]

    if a0 is not None:
        a0 = np.asarray(a0)
        if a0.shape != (n, n):
            raise ValueError(f"a0 must have shape {(n,n)}, got {a0.shape}")

    z = float(x) + 1j * float(eps)

    if a0 is None:
        G = cauchy_matrix_semicircle(z, A, G0=G0, tol=tol, maxiter=maxiter)
    else:
        G = cauchy_biased_matrix_semicircle(z, a0, A, G0=G0, tol=tol, maxiter=maxiter)

    m = np.trace(G) / n
    f = (-1.0 / np.pi) * np.imag(m)
    return float(f)

biased_matrix_semicircle_density(x, a0, A, eps=0.01, G0=None, tol=1e-10, maxiter=10000)

Scalar density of the biased matrix semicircle \(S = a_0 + \sum_i A_i \otimes X_i\).

Convenience wrapper: $$ f_{a_0}(x) \;=\; -\frac{1}{\pi}\,\Im!\left(\frac{1}{n}\, \mathrm{tr}\,G_{a_0+X}(x+i\varepsilon)\right). $$

Calls matrix_semicircle_density(x, A, eps, G0, tol, maxiter, a0=a0).

Source code in src/free_matrix_laws/transforms.py

def biased_matrix_semicircle_density(
    x: float,
    a0,
    A,
    eps: float = 1e-2,
    G0=None,
    tol: float = 1e-10,
    maxiter: int = 10_000,
) -> float:
    r'''
    Scalar density of the biased matrix semicircle
    $S = a_0 + \sum_i A_i \otimes X_i$.

    Convenience wrapper:
    $$ f_{a_0}(x) \;=\; -\frac{1}{\pi}\,\Im\!\left(\frac{1}{n}\,
       \mathrm{tr}\,G_{a_0+X}(x+i\varepsilon)\right). $$

    Calls ``matrix_semicircle_density(x, A, eps, G0, tol, maxiter, a0=a0)``.
    '''
    return matrix_semicircle_density(x, A, eps=eps, G0=G0, tol=tol, maxiter=maxiter, a0=a0)

polynomial_density(x, a0, A, *, eps=0.01, eps_reg=None, block_size=1, G0=None, tol=1e-10, maxiter=10000, return_info=False)

Scalar density of a self-adjoint polynomial \(p(X_1,\dots,X_s)\) of free semicircular variables, via self-adjoint linearization.

Evaluates at \(z=x+i\,\varepsilon\) and computes the regularized fixed point \(G(z,b_\varepsilon(z))\) via :func:cauchy_polynomial.

The scalar Cauchy transform is extracted from the distinguished corner: $$ m_p(z) \approx \frac{1}{k}\,\mathrm{tr}\, \big(G(z,b_\varepsilon(z))\big)_{11}, $$ and the density is approximated by $$ f(x) \approx -\frac{1}{\pi}\,\Im\, m_p(x+i\varepsilon). $$

Parameters:

Name	Type	Description	Default
x	float	Real evaluation point.	required
a0	(n,n) array	Constant term of the self-adjoint linearization.	required
A	sequence of (n,n) arrays or stacked array (s,n,n)	Coefficients defining \(\eta(B)=\sum_i A_i B A_i^\ast\).	required
eps	float	Imaginary offset in \(z=x+i\,\varepsilon\).	1e-2
eps_reg	float	Regularization in \(\Lambda_\varepsilon(z)\). If None, uses eps.	None
block_size	int	Size \(k\) of the distinguished top-left block.	1
G0	(n,n) array	Warm start for the solver.	None
tol	float	Relative tolerance.	1e-10
maxiter	int	Maximum iterations.	10000
return_info	bool	If True, also return solver diagnostics.	False

Returns:

Type	Description
float	Approximation to the density \(f(x)\).

Source code in src/free_matrix_laws/transforms.py

def polynomial_density(
    x: float,
    a0: np.ndarray,
    A,
    *,
    eps: float = 1e-2,
    eps_reg: float | None = None,
    block_size: int = 1,
    G0: np.ndarray | None = None,
    tol: float = 1e-10,
    maxiter: int = 10_000,
    return_info: bool = False,
) -> float:
    r'''
    Scalar density of a self-adjoint polynomial $p(X_1,\dots,X_s)$ of free
    semicircular variables, via self-adjoint linearization.

    Evaluates at $z=x+i\,\varepsilon$ and computes the regularized fixed point
    $G(z,b_\varepsilon(z))$ via :func:`cauchy_polynomial`.

    The scalar Cauchy transform is extracted from the distinguished corner:
    $$
      m_p(z) \approx \frac{1}{k}\,\mathrm{tr}\,
          \big(G(z,b_\varepsilon(z))\big)_{11},
    $$
    and the density is approximated by
    $$
      f(x) \approx -\frac{1}{\pi}\,\Im\, m_p(x+i\varepsilon).
    $$

    Parameters
    ----------
    x : float
        Real evaluation point.
    a0 : (n,n) array
        Constant term of the self-adjoint linearization.
    A : sequence of (n,n) arrays or stacked array (s,n,n)
        Coefficients defining $\eta(B)=\sum_i A_i B A_i^\ast$.
    eps : float, default 1e-2
        Imaginary offset in $z=x+i\,\varepsilon$.
    eps_reg : float, optional
        Regularization in $\Lambda_\varepsilon(z)$. If None, uses eps.
    block_size : int, default 1
        Size $k$ of the distinguished top-left block.
    G0 : (n,n) array, optional
        Warm start for the solver.
    tol : float, default 1e-10
        Relative tolerance.
    maxiter : int, default 10000
        Maximum iterations.
    return_info : bool, default False
        If True, also return solver diagnostics.

    Returns
    -------
    float
        Approximation to the density $f(x)$.
    '''
    if eps <= 0:
        raise ValueError("eps must be > 0.")
    z = float(x) + 1j * float(eps)
    if eps_reg is None:
        eps_reg = float(eps)

    G, info = cauchy_polynomial(
        z, a0, A,
        eps_reg=eps_reg, block_size=block_size,
        G0=G0, tol=tol, maxiter=maxiter,
        return_info=True,
    )

    k = int(block_size)
    G11 = G[:k, :k]
    m = np.trace(G11) / k
    f = (-1.0 / np.pi) * np.imag(m)

    if return_info:
        info = dict(info)
        info["z"] = z
        info["m"] = complex(m)
        info["density"] = float(f)
        return float(f), info

    return float(f)

semicircle_density_scalar(x, c=1.0)

Classical (scalar) Wigner semicircle density with variance \(c>0\).

\[ f(x) \;=\; \frac{1}{2\pi c}\,\sqrt{\,4c - x^2\,}\, \mathbf 1_{\{|x|\le 2\sqrt{c}\}}. \]

Parameters:

Name	Type	Description	Default
x	float or array_like		required
c	float	Variance parameter (\(c>0\), so radius is \(2\sqrt{c}\)).	1.0

Returns:

Type	Description
float or ndarray

Source code in src/free_matrix_laws/transforms.py

def semicircle_density_scalar(x, c: float = 1.0):
    r'''
    Classical (scalar) Wigner semicircle density with variance $c>0$.

    $$
      f(x) \;=\; \frac{1}{2\pi c}\,\sqrt{\,4c - x^2\,}\,
      \mathbf 1_{\{|x|\le 2\sqrt{c}\}}.
    $$

    Parameters
    ----------
    x : float or array_like
    c : float, default 1.0
        Variance parameter ($c>0$, so radius is $2\sqrt{c}$).

    Returns
    -------
    float or ndarray
    '''
    if c <= 0:
        raise ValueError("c must be > 0")
    x_arr = np.asarray(x, dtype=float)
    inside = 4.0 * c - x_arr**2
    y = np.where(inside > 0.0, (1.0 / (2.0 * np.pi * c)) * np.sqrt(inside), 0.0)
    return y if x_arr.ndim else float(y)

semicircle_cauchy_scalar(z, c=1.0)

Scalar Cauchy (Stieltjes) transform of the Wigner semicircle law with variance \(c>0\).

\[G(z) = \frac{z - \sqrt{z^2 - 4c}}{2c}\]

The square-root branch is chosen so that \(\Im z>0 \Rightarrow \Im G(z)<0\).

Source code in src/free_matrix_laws/transforms.py

def semicircle_cauchy_scalar(z, c: float = 1.0):
    r"""
    Scalar Cauchy (Stieltjes) transform of the Wigner semicircle law
    with variance $c>0$.

    $$G(z) = \frac{z - \sqrt{z^2 - 4c}}{2c}$$

    The square-root branch is chosen so that $\Im z>0 \Rightarrow \Im G(z)<0$.
    """
    if c <= 0:
        raise ValueError("c must be > 0")

    z_arr = np.asarray(z, dtype=np.complex128)
    disc = np.sqrt(z_arr**2 - 4.0 * c)
    disc = np.where(disc.imag * z_arr.imag < 0, -disc, disc)

    G = (z_arr - disc) / (2.0 * c)
    return G if z_arr.ndim else complex(G)

cauchy_kronecker(w, a, cauchy_scalar=None, eps=1e-08)

Cauchy transform of a Kronecker-product random variable \(a \otimes X\).

Computes $$ G_a(w) \;=\; E!\big[(w - a \otimes X)^{-1}\big], $$ where \(X\) is a scalar random variable with known Cauchy (Stieltjes) transform \(G_X(z)\), and \(a\) is an \(n \times n\) deterministic matrix.

Method. Regularize \(a\) to \(\tilde a = a + i\varepsilon I\) so that it is invertible, then diagonalize \(\tilde a^{-1} w = V\,\mathrm{diag}(\mu_1,\dots,\mu_n)\,V^{-1}\). The expectation reduces to $$ G_a(w) \;=\; V\,\mathrm{diag}!\big(G_X(\mu_1),\dots,G_X(\mu_n)\big)\, V^{-1}\,\tilde a^{-1}, $$ where \(G_X\) is the scalar Cauchy transform passed as cauchy_scalar.

This is \(O(n^3)\) (one eigendecomposition + two solves) regardless of the underlying distribution of \(X\).

Parameters:

Name	Type	Description	Default
w	(n, n) ndarray	Spectral parameter matrix (should have \(\Im w \ne 0\) in some sense, e.g. \(w = (x + i\varepsilon)\,I\)).	required
a	(n, n) ndarray	The deterministic matrix in \(a \otimes X\).	required
cauchy_scalar	callable	Scalar Cauchy transform \(G_X(z)\), accepting a complex array and returning a complex array of the same shape. Default: :func:semicircle_cauchy_scalar (standard Wigner law).	None
eps	float	Regularization: replaces \(a\) by \(a + i\varepsilon I\) to handle singular or near-singular \(a\).	1e-8

Returns:

Type	Description
(n, n) ndarray (complex)	The matrix-valued Cauchy transform \(G_a(w)\).

See Also

cauchy_kronecker_semicircle : Convenience alias with cauchy_scalar=semicircle_cauchy_scalar. h_kronecker : Subordination h-function for the same model. cauchy_matrix_semicircle : General case \(\sum_i A_i \otimes X_i\) (fixed-point iteration).

Examples:

>>> import numpy as np
>>> from free_matrix_laws import cauchy_kronecker, semicircle_cauchy_scalar
>>> w = (0.5 + 0.01j) * np.eye(2)
>>> a = np.array([[1.0, 0.2], [0.2, 0.8]])
>>> G = cauchy_kronecker(w, a)  # semicircle by default
>>> G.shape
(2, 2)

Use a custom scalar Cauchy transform (e.g. Marchenko–Pastur):

>>> def cauchy_mp(z, gamma=1.0):
...     # Marchenko-Pastur Cauchy transform
...     disc = np.sqrt((z - (1+np.sqrt(gamma))**2) *
...                    (z - (1-np.sqrt(gamma))**2))
...     disc = np.where(disc.imag * z.imag < 0, -disc, disc)
...     return (z - gamma + 1 - disc) / (2 * gamma * z)
>>> G_mp = cauchy_kronecker(w, a, cauchy_scalar=cauchy_mp)

Source code in src/free_matrix_laws/transforms.py

def cauchy_kronecker(
    w: np.ndarray,
    a: np.ndarray,
    cauchy_scalar: Callable = None,
    eps: float = 1e-8,
) -> np.ndarray:
    r'''
    Cauchy transform of a Kronecker-product random variable $a \otimes X$.

    Computes
    $$
      G_a(w) \;=\; E\!\big[(w - a \otimes X)^{-1}\big],
    $$
    where $X$ is a scalar random variable with known Cauchy (Stieltjes) transform
    $G_X(z)$, and $a$ is an $n \times n$ deterministic matrix.

    **Method.** Regularize $a$ to $\tilde a = a + i\varepsilon I$ so that it is
    invertible, then diagonalize
    $\tilde a^{-1} w = V\,\mathrm{diag}(\mu_1,\dots,\mu_n)\,V^{-1}$.
    The expectation reduces to
    $$
      G_a(w) \;=\; V\,\mathrm{diag}\!\big(G_X(\mu_1),\dots,G_X(\mu_n)\big)\,
                    V^{-1}\,\tilde a^{-1},
    $$
    where $G_X$ is the scalar Cauchy transform passed as ``cauchy_scalar``.

    This is $O(n^3)$ (one eigendecomposition + two solves) regardless of
    the underlying distribution of $X$.

    Parameters
    ----------
    w : (n, n) ndarray
        Spectral parameter matrix (should have $\Im w \ne 0$ in some sense,
        e.g. $w = (x + i\varepsilon)\,I$).
    a : (n, n) ndarray
        The deterministic matrix in $a \otimes X$.
    cauchy_scalar : callable, optional
        Scalar Cauchy transform $G_X(z)$, accepting a complex array and
        returning a complex array of the same shape.
        Default: :func:`semicircle_cauchy_scalar` (standard Wigner law).
    eps : float, default 1e-8
        Regularization: replaces $a$ by $a + i\varepsilon I$ to handle
        singular or near-singular $a$.

    Returns
    -------
    (n, n) ndarray (complex)
        The matrix-valued Cauchy transform $G_a(w)$.

    See Also
    --------
    cauchy_kronecker_semicircle :
        Convenience alias with ``cauchy_scalar=semicircle_cauchy_scalar``.
    h_kronecker :
        Subordination h-function for the same model.
    cauchy_matrix_semicircle :
        General case $\sum_i A_i \otimes X_i$ (fixed-point iteration).

    Examples
    --------
    >>> import numpy as np
    >>> from free_matrix_laws import cauchy_kronecker, semicircle_cauchy_scalar
    >>> w = (0.5 + 0.01j) * np.eye(2)
    >>> a = np.array([[1.0, 0.2], [0.2, 0.8]])
    >>> G = cauchy_kronecker(w, a)  # semicircle by default
    >>> G.shape
    (2, 2)

    Use a custom scalar Cauchy transform (e.g. Marchenko–Pastur):

    >>> def cauchy_mp(z, gamma=1.0):
    ...     # Marchenko-Pastur Cauchy transform
    ...     disc = np.sqrt((z - (1+np.sqrt(gamma))**2) *
    ...                    (z - (1-np.sqrt(gamma))**2))
    ...     disc = np.where(disc.imag * z.imag < 0, -disc, disc)
    ...     return (z - gamma + 1 - disc) / (2 * gamma * z)
    >>> G_mp = cauchy_kronecker(w, a, cauchy_scalar=cauchy_mp)
    '''
    if cauchy_scalar is None:
        cauchy_scalar = semicircle_cauchy_scalar

    w = np.asarray(w, dtype=complex)
    a = np.asarray(a, dtype=complex)
    if w.ndim != 2 or w.shape[0] != w.shape[1]:
        raise ValueError(f"w must be square; got {w.shape!r}")
    if a.shape != w.shape:
        raise ValueError(f"a must have same shape as w {w.shape!r}; got {a.shape!r}")
    if eps < 0:
        raise ValueError("eps must be >= 0")

    n = w.shape[0]
    a_reg = a + 1j * eps * np.eye(n, dtype=complex)
    a_reg_inv = la.inv(a_reg)

    mu, V = la.eig(a_reg_inv @ w)
    G_mu = cauchy_scalar(mu)  # vectorized over eigenvalues

    return V @ np.diag(G_mu) @ la.inv(V) @ a_reg_inv

h_kronecker(w, a, cauchy_scalar=None, eps=1e-08)

Subordination h-function for a Kronecker-product random variable \(a \otimes X\).

$$ h_a(w) \;=\; G_a(w)^{-1} \;-\; w, $$ where \(G_a(w)\) is computed by :func:cauchy_kronecker.

Parameters:

Name	Type	Description	Default
w	(n, n) ndarray	Spectral parameter matrix.	required
a	(n, n) ndarray	Deterministic matrix in \(a \otimes X\).	required
cauchy_scalar	callable	Scalar Cauchy transform \(G_X(z)\). Default: :func:semicircle_cauchy_scalar.	None
eps	float	Regularization parameter.	1e-8

Returns:

Type	Description
(n, n) ndarray (complex)

See Also

h_kronecker_semicircle : Convenience alias for the semicircle case. cauchy_kronecker : The underlying Cauchy transform.

Source code in src/free_matrix_laws/transforms.py

def h_kronecker(
    w: np.ndarray,
    a: np.ndarray,
    cauchy_scalar: Callable = None,
    eps: float = 1e-8,
) -> np.ndarray:
    r'''
    Subordination h-function for a Kronecker-product random variable $a \otimes X$.

    $$
      h_a(w) \;=\; G_a(w)^{-1} \;-\; w,
    $$
    where $G_a(w)$ is computed by :func:`cauchy_kronecker`.

    Parameters
    ----------
    w : (n, n) ndarray
        Spectral parameter matrix.
    a : (n, n) ndarray
        Deterministic matrix in $a \otimes X$.
    cauchy_scalar : callable, optional
        Scalar Cauchy transform $G_X(z)$.
        Default: :func:`semicircle_cauchy_scalar`.
    eps : float, default 1e-8
        Regularization parameter.

    Returns
    -------
    (n, n) ndarray (complex)

    See Also
    --------
    h_kronecker_semicircle :
        Convenience alias for the semicircle case.
    cauchy_kronecker :
        The underlying Cauchy transform.
    '''
    G = cauchy_kronecker(w, a, cauchy_scalar=cauchy_scalar, eps=eps)
    return la.inv(G) - w

cauchy_kronecker_semicircle(w, a, eps=1e-08)

Cauchy transform of the Kronecker-product semicircle \(a \otimes X\).

Convenience alias for cauchy_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps).

See :func:cauchy_kronecker for full documentation.

Parameters:

Name	Type	Description	Default
w	(n, n) ndarray	Spectral parameter matrix.	required
a	(n, n) ndarray	Deterministic matrix in \(a \otimes X\).	required
eps	float	Regularization parameter.	1e-8

Returns:

Type	Description
(n, n) ndarray (complex)

Source code in src/free_matrix_laws/transforms.py

def cauchy_kronecker_semicircle(
    w: np.ndarray,
    a: np.ndarray,
    eps: float = 1e-8,
) -> np.ndarray:
    r'''
    Cauchy transform of the Kronecker-product semicircle $a \otimes X$.

    Convenience alias for
    ``cauchy_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps)``.

    See :func:`cauchy_kronecker` for full documentation.

    Parameters
    ----------
    w : (n, n) ndarray
        Spectral parameter matrix.
    a : (n, n) ndarray
        Deterministic matrix in $a \otimes X$.
    eps : float, default 1e-8
        Regularization parameter.

    Returns
    -------
    (n, n) ndarray (complex)
    '''
    return cauchy_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps)

h_kronecker_semicircle(w, a, eps=1e-08)

Subordination h-function for the Kronecker-product semicircle \(a \otimes X\).

Convenience alias for h_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps).

See :func:h_kronecker for full documentation.

Parameters:

Name	Type	Description	Default
w	(n, n) ndarray	Spectral parameter matrix.	required
a	(n, n) ndarray	Deterministic matrix in \(a \otimes X\).	required
eps	float	Regularization parameter.	1e-8

Returns:

Type	Description
(n, n) ndarray (complex)

Source code in src/free_matrix_laws/transforms.py

def h_kronecker_semicircle(
    w: np.ndarray,
    a: np.ndarray,
    eps: float = 1e-8,
) -> np.ndarray:
    r'''
    Subordination h-function for the Kronecker-product semicircle $a \otimes X$.

    Convenience alias for
    ``h_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps)``.

    See :func:`h_kronecker` for full documentation.

    Parameters
    ----------
    w : (n, n) ndarray
        Spectral parameter matrix.
    a : (n, n) ndarray
        Deterministic matrix in $a \otimes X$.
    eps : float, default 1e-8
        Regularization parameter.

    Returns
    -------
    (n, n) ndarray (complex)
    '''
    return h_kronecker(w, a, cauchy_scalar=semicircle_cauchy_scalar, eps=eps)

subordination_kronecker(b, a1, a2, cauchy_scalar_x=None, cauchy_scalar_y=None, eps=0.0001, tol=1e-08, maxiter=10000, return_info=False)

Subordination function \(\omega_1(b)\) for the free additive convolution of two Kronecker-product random variables \(a_1 \otimes X\) and \(a_2 \otimes Y\).

Computes the fixed point of the map $$ w \;\mapsto\; h_Y!\big(h_X(w) + b\big) + b, $$ where \(h_X(w) = G_{a_1 \otimes X}(w)^{-1} - w\) and \(h_Y(w) = G_{a_2 \otimes Y}(w)^{-1} - w\) are the subordination h-functions computed via :func:h_kronecker.

After convergence, the Cauchy transform of the sum \(p(X,Y)\) (as encoded by the linearization \(b = \Lambda_\varepsilon(z) - a_0\)) can be recovered as $$ G_{X+Y}(b) \;=\; G_{a_1 \otimes X}!\big(\omega_1(b)\big). $$

Parameters:

Name	Type	Description	Default
b	(n, n) ndarray	The "driving matrix," typically \(b = \Lambda_\varepsilon(z) - a_0\) for a linearization \(L_p = a_0 + a_1 \otimes X + a_2 \otimes Y\).	required
a1	(n, n) ndarray	Deterministic matrix for the first variable (\(a_1 \otimes X\)).	required
a2	(n, n) ndarray	Deterministic matrix for the second variable (\(a_2 \otimes Y\)).	required
cauchy_scalar_x	callable	Scalar Cauchy transform \(G_X(z)\) for the first variable. Default: :func:semicircle_cauchy_scalar.	None
cauchy_scalar_y	callable	Scalar Cauchy transform \(G_Y(z)\) for the second variable. Default: :func:semicircle_cauchy_scalar.	None
eps	float	Regularization for the Kronecker h-functions. Larger than the default in :func:cauchy_kronecker because linearization matrices \(a_1, a_2\) are typically rank-deficient.	1e-4
tol	float	Convergence tolerance (Frobenius norm of \(w_{k+1} - w_k\)). Note: with eps=1e-4, achievable accuracy is roughly \(O(\varepsilon)\), so tighter tolerances may not be reached.	1e-8
maxiter	int	Maximum iterations.	10000
return_info	bool	If True, also return a dict with diagnostics.	False

Returns:

Name	Type	Description
omega	(n, n) ndarray (complex)	The subordination function \(\omega_1(b)\).
info	dict (only if return_info=True)	Keys: iters, last_diff.

See Also

h_kronecker : The h-function used at each step. cauchy_kronecker : To recover \(G_{X+Y}(b)\) from \(\omega_1(b)\).

Examples:

Anticommutator of two free semicircles via subordination:

>>> import numpy as np
>>> from free_matrix_laws import (
...     subordination_kronecker, cauchy_kronecker_semicircle, lambda_eps
... )
>>> A0 = np.array([[0, 0, 0], [0, 0, -1], [0, -1, 0]])
>>> A1 = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
>>> A2 = np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]])
>>> z = 0.5 + 0.01j
>>> b = lambda_eps(z, 3) - A0
>>> omega = subordination_kronecker(b, A1, A2)
>>> G = cauchy_kronecker_semicircle(omega, A1)
>>> density = -G[0, 0].imag / np.pi

Anticommutator with different distributions (e.g. free Poisson):

>>> def cauchy_poisson(z, lam=4.0):
...     a = (1 - np.sqrt(lam))**2
...     b = (1 + np.sqrt(lam))**2
...     disc = np.sqrt((z - a) * (z - b))
...     disc = np.where(disc.imag * z.imag < 0, -disc, disc)
...     return (1 + z - lam - disc) / (2 * z)
>>> omega = subordination_kronecker(b, A1, A2,
...     cauchy_scalar_x=cauchy_poisson, cauchy_scalar_y=cauchy_poisson)

References

• S. Belinschi, T. Mai, R. Speicher, Analytic subordination theory of operator-valued free additive convolution and the solution of a general random matrix problem, J. reine angew. Math. 732 (2017), 21–53.

Source code in src/free_matrix_laws/transforms.py

def subordination_kronecker(
    b: np.ndarray,
    a1: np.ndarray,
    a2: np.ndarray,
    cauchy_scalar_x: Callable = None,
    cauchy_scalar_y: Callable = None,
    eps: float = 1e-4,
    tol: float = 1e-8,
    maxiter: int = 10_000,
    return_info: bool = False,
) -> np.ndarray:
    r'''
    Subordination function $\omega_1(b)$ for the free additive convolution
    of two Kronecker-product random variables $a_1 \otimes X$ and $a_2 \otimes Y$.

    Computes the fixed point of the map
    $$
      w \;\mapsto\; h_Y\!\big(h_X(w) + b\big) + b,
    $$
    where $h_X(w) = G_{a_1 \otimes X}(w)^{-1} - w$ and
    $h_Y(w) = G_{a_2 \otimes Y}(w)^{-1} - w$ are the subordination h-functions
    computed via :func:`h_kronecker`.

    After convergence, the Cauchy transform of the sum $p(X,Y)$ (as encoded
    by the linearization $b = \Lambda_\varepsilon(z) - a_0$) can be recovered as
    $$
      G_{X+Y}(b) \;=\; G_{a_1 \otimes X}\!\big(\omega_1(b)\big).
    $$

    Parameters
    ----------
    b : (n, n) ndarray
        The "driving matrix," typically $b = \Lambda_\varepsilon(z) - a_0$
        for a linearization $L_p = a_0 + a_1 \otimes X + a_2 \otimes Y$.
    a1 : (n, n) ndarray
        Deterministic matrix for the first variable ($a_1 \otimes X$).
    a2 : (n, n) ndarray
        Deterministic matrix for the second variable ($a_2 \otimes Y$).
    cauchy_scalar_x : callable, optional
        Scalar Cauchy transform $G_X(z)$ for the first variable.
        Default: :func:`semicircle_cauchy_scalar`.
    cauchy_scalar_y : callable, optional
        Scalar Cauchy transform $G_Y(z)$ for the second variable.
        Default: :func:`semicircle_cauchy_scalar`.
    eps : float, default 1e-4
        Regularization for the Kronecker h-functions. Larger than the
        default in :func:`cauchy_kronecker` because linearization matrices
        $a_1, a_2$ are typically rank-deficient.
    tol : float, default 1e-8
        Convergence tolerance (Frobenius norm of $w_{k+1} - w_k$).
        Note: with ``eps=1e-4``, achievable accuracy is roughly $O(\varepsilon)$,
        so tighter tolerances may not be reached.
    maxiter : int, default 10000
        Maximum iterations.
    return_info : bool, default False
        If True, also return a dict with diagnostics.

    Returns
    -------
    omega : (n, n) ndarray (complex)
        The subordination function $\omega_1(b)$.
    info : dict (only if return_info=True)
        Keys: ``iters``, ``last_diff``.

    See Also
    --------
    h_kronecker : The h-function used at each step.
    cauchy_kronecker : To recover $G_{X+Y}(b)$ from $\omega_1(b)$.

    Examples
    --------
    Anticommutator of two free semicircles via subordination:

    >>> import numpy as np
    >>> from free_matrix_laws import (
    ...     subordination_kronecker, cauchy_kronecker_semicircle, lambda_eps
    ... )
    >>> A0 = np.array([[0, 0, 0], [0, 0, -1], [0, -1, 0]])
    >>> A1 = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
    >>> A2 = np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]])
    >>> z = 0.5 + 0.01j
    >>> b = lambda_eps(z, 3) - A0
    >>> omega = subordination_kronecker(b, A1, A2)
    >>> G = cauchy_kronecker_semicircle(omega, A1)
    >>> density = -G[0, 0].imag / np.pi

    Anticommutator with different distributions (e.g. free Poisson):

    >>> def cauchy_poisson(z, lam=4.0):
    ...     a = (1 - np.sqrt(lam))**2
    ...     b = (1 + np.sqrt(lam))**2
    ...     disc = np.sqrt((z - a) * (z - b))
    ...     disc = np.where(disc.imag * z.imag < 0, -disc, disc)
    ...     return (1 + z - lam - disc) / (2 * z)
    >>> omega = subordination_kronecker(b, A1, A2,
    ...     cauchy_scalar_x=cauchy_poisson, cauchy_scalar_y=cauchy_poisson)

    References
    ----------
    * S. Belinschi, T. Mai, R. Speicher, *Analytic subordination theory of
      operator-valued free additive convolution and the solution of a general
      random matrix problem*, J. reine angew. Math. **732** (2017), 21–53.
    '''
    if cauchy_scalar_x is None:
        cauchy_scalar_x = semicircle_cauchy_scalar
    if cauchy_scalar_y is None:
        cauchy_scalar_y = semicircle_cauchy_scalar

    b = np.asarray(b, dtype=complex)
    a1 = np.asarray(a1, dtype=complex)
    a2 = np.asarray(a2, dtype=complex)
    n = b.shape[0]

    W = 1j * np.eye(n, dtype=complex)  # initialization

    last_diff = np.inf
    for k in range(1, maxiter + 1):
        W1 = h_kronecker(W, a1, cauchy_scalar=cauchy_scalar_x, eps=eps) + b
        W_new = h_kronecker(W1, a2, cauchy_scalar=cauchy_scalar_y, eps=eps) + b

        diff = la.norm(W_new - W, 'fro')
        last_diff = diff

        if diff <= tol:
            W = W_new
            break
        W = W_new

    if return_info:
        return W, {"iters": k, "last_diff": float(last_diff)}
    return W

lambda_eps(z, n, eps=1e-06, block_size=1)

Regularized spectral parameter for polynomial linearizations.

\[ \Lambda_\varepsilon(z) = \begin{bmatrix} z\, I_{k} & 0 \\ 0 & i\varepsilon\, I_{n-k} \end{bmatrix}, \qquad k = \texttt{block\_size}. \]

Parameters:

Name	Type	Description	Default
z	complex	Spectral parameter.	required
n	int	Matrix size.	required
eps	float	Regularization \(\varepsilon > 0\) on the lower block.	``1e-6``
block_size	int	Size \(k\) of the top-left block carrying \(z\).	``1``

Returns:

Type	Description
(n, n) ndarray, complex

Source code in src/free_matrix_laws/transforms.py

def lambda_eps(z: complex, n: int, eps: float = 1e-6, block_size: int = 1) -> np.ndarray:
    r'''
    Regularized spectral parameter for polynomial linearizations.

    $$
      \Lambda_\varepsilon(z)
      =
      \begin{bmatrix}
        z\, I_{k} & 0 \\
        0 & i\varepsilon\, I_{n-k}
      \end{bmatrix},
    \qquad k = \texttt{block\_size}.
    $$

    Parameters
    ----------
    z : complex
        Spectral parameter.
    n : int
        Matrix size.
    eps : float, default ``1e-6``
        Regularization $\varepsilon > 0$ on the lower block.
    block_size : int, default ``1``
        Size $k$ of the top-left block carrying $z$.

    Returns
    -------
    (n, n) ndarray, complex
    '''
    return _lambda_eps(z, n, eps_reg=eps, block_size=block_size)