Skip to content

residue_id ¤

Stage 2: Residue identification via weighted least squares.

Ported from vectfit4.m lines 542–750 (SINTEF MFT-NNLS, B. Gustavsen, 2026).

The core computation is a per-element LS solve: for each of the nnn upper-triangle elements, find residues C, and optionally D and E terms, given fixed poles. This is the primary vmap target: each element is independent.

Functions:

Name Description
identify_residues

Identify residues (and D, E terms) for all upper-triangle elements.

identify_residues ¤

identify_residues(
    f: ndarray, s: ndarray, poles: ndarray, weight: ndarray, opts
) -> tuple[ndarray, ndarray, ndarray]

Identify residues (and D, E terms) for all upper-triangle elements.

Returns:

Name Type Description
C ndarray

(nnn, N) complex — residues for each element and each pole.
D_vec ndarray

(nnn,) float — constant terms (zero if opts.asymp == 1).
E_vec ndarray

(nnn,) float — linear terms (zero if opts.asymp < 3).
Source code in circulax/fitting/residue_id.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def identify_residues(
    f: jnp.ndarray,  # (nnn, Ns) complex — upper-triangle elements
    s: jnp.ndarray,  # (Ns,) complex
    poles: np.ndarray,  # (N,) complex NumPy — current (fixed) poles
    weight: jnp.ndarray,  # (nnn, Ns) or (1, Ns) — fitting weights
    opts,  # FitOptions
) -> tuple[jnp.ndarray, jnp.ndarray, jnp.ndarray]:
    """Identify residues (and D, E terms) for all upper-triangle elements.

    Returns:
        C:     (nnn, N) complex — residues for each element and each pole.
        D_vec: (nnn,) float    — constant terms (zero if opts.asymp == 1).
        E_vec: (nnn,) float    — linear terms   (zero if opts.asymp < 3).

    """
    N = len(poles)
    Ns = s.shape[0]
    nnn = f.shape[0]
    offs = opts.asymp - 1  # number of extra columns (0, 1, or 2)

    cindex = compute_cindex(poles)
    Dk_uw = build_dk(s, jnp.array(poles), cindex, offs)  # (Ns, N+offs)

    # Broadcast weight to (nnn, Ns)
    w = jnp.broadcast_to(weight, (nnn, Ns)).real  # weights are real

    # vmap the per-element LS solve over the nnn elements
    _solve_vmap = jax.vmap(
        lambda f_n, w_n: _residue_ls_one(f_n, w_n, Dk_uw, s, offs),
        in_axes=(0, 0),
    )
    x_all = _solve_vmap(f, w)  # (nnn, N+offs) real

    # Split into residue coefficients and D, E
    C_real = np.array(x_all[:, :N])  # (nnn, N) — real representation

    # Convert each element's residues from real to complex
    C_complex = np.stack(
        [_real_to_complex_residues(C_real[n], cindex, N) for n in range(nnn)],
        axis=0,
    )  # (nnn, N) complex

    # Extract D and E from the extra columns
    if opts.asymp >= 2:
        D_vec = jnp.real(x_all[:, N])  # (nnn,) — D term
    else:
        D_vec = jnp.zeros(nnn)

    if opts.asymp >= 3:
        E_vec = jnp.real(x_all[:, N + 1])  # (nnn,) — E term
    else:
        E_vec = jnp.zeros(nnn)

    return jnp.array(C_complex), D_vec, E_vec