An autonomous CPU implementation of the Janashia-Lagvilava matrix spectral factorization algorithm using NumPy and SciPy.
Translation and status notice: This repository is an independent Python translation of the original JLAlgorithm MATLAB repository. The translation was produced with the assistance of a large language model and should be treated as a modified work based on the original project. It has not been formally validated against every original implementation path. No performance, accuracy, numerical stability, or fitness-for-purpose guarantees are provided. Always inspect the reported diagnostics and validate results for your own data and application.
The original MATLAB repository and this Python translation are released under the GNU General Public License, version 3.0. See LICENSE for the complete license text. This translation does not change the licensing or patent caveats associated with the underlying algorithm.
The package accepts matrix-polynomial coefficients with shape (d, d, degree + 1) and returns causal factor coefficients in the same layout. The implementation is independent of MATLAB at runtime and in its test suite.
import numpy as np
from jl_factorization import block_factorization
coefficients = np.zeros((4, 4, 3), dtype=float)
coefficients[:, :, 0] = np.eye(4)
coefficients[:, :, 1:] = 0.02 * np.random.default_rng(0).standard_normal((4, 4, 2))
result = block_factorization(coefficients, fft_power=5)
factor_coefficients = result.coefficients
print(result.diagnostics.frequency_error)block_factorization currently requires a power-of-two matrix dimension and an FFT size large enough for the product polynomial. scalar_factorization is also available for the original iterative path. Both return a FactorizationResult containing the coefficient tensor and diagnostics.
The hierarchical block path is the validated implementation target and currently reconstructs small and medium deterministic spectra to machine precision. The scalar path is present as a direct reference implementation, but its transform parity is still being refined; callers should inspect result.diagnostics.frequency_error rather than assuming the scalar path has the same accuracy yet.
python -m pip install -e ".[test]"
python -m pytestRun the complete example from this project directory:
python examples/demo_factorization.pyTo use the package from another project, install it by referring to this directory, then import the public solver:
python -m pip install -e /path/to/JLAlgorithm-pythonfrom jl_factorization import block_factorization
result = block_factorization(coefficients, fft_power=5)
factor_coefficients = result.coefficientsThe initial release targets NumPy/SciPy CPU execution. GPU and parallel backends are intentionally separate future extensions.
The underlying algorithm is subject to the licensing and patent caveats described by the original project. Users are responsible for evaluating any patent or commercial-use obligations.