Add support for deflation in arnoldi_decomposition #14
Conversation
We also add a test to make sure it works in the special case of early break (== finding an invariant subspace) at the first iteration. This will often happen when doing restarts based methods.
This opens the door to deflation/locking.
We ensure that arnoldi invariant are kept when initializing from a set of eigen-pairs.
There was a problem hiding this comment.
Pull Request Overview
This PR adds support for deflation in the arnoldi_decomposition function, enabling it to handle "lucky breaks" where the Krylov basis converges earlier than expected. The implementation includes explicit restart capabilities and improved parameter handling.
- Adds deflation logic to handle early convergence scenarios ("lucky breaks")
- Implements explicit restart functionality with
naive_explicit_restarts - Introduces flexible start/stop capabilities for the decomposition process
Reviewed Changes
Copilot reviewed 6 out of 6 changed files in this pull request and generated 5 comments.
Show a summary per file
| File | Description |
|---|---|
src/arnoldi/decomposition.py |
Adds deflation support, start_dim parameter, and default tolerance handling |
src/arnoldi/explicit_restarts.py |
New module implementing naive explicit restart algorithm |
tests/test_decomposition.py |
Comprehensive tests for new deflation and restart features |
tests/test_explicit_restarts.py |
Test suite for explicit restart functionality |
tests/common.py |
Shared test constants for retry configuration |
Makefile |
Development tooling for linting and testing |
Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.
| """Create the basis vector e_k in R^n, aka e_k is (n,), and with e_k[k-1] = | ||
| 1 | ||
| """ |
There was a problem hiding this comment.
The updated comment is inconsistent with the function implementation. The function sets ret[k-1] = 1, but the comment suggests this creates e_k[k-1] = 1, which is confusing since e_k typically refers to the k-th standard basis vector where the k-th element (1-indexed) is 1.
| """Create the basis vector e_k in R^n, aka e_k is (n,), and with e_k[k-1] = | |
| 1 | |
| """ | |
| """ | |
| Create the k-th standard basis vector in R^n (of shape (n,)), where the entry at index (k-1) | |
| is 1 and all other entries are 0. (Note: k is 1-based, as in mathematical notation; Python uses 0-based indexing.) | |
| """ |
| # the arnoldi basis V is orthonormal | ||
| np.testing.assert_allclose( | ||
| V.conj().T @ V, np.eye(m + 1), rtol=RTOL, atol=ATOL | ||
| V_m.conj().T @ V_m, np.eye(m), rtol=RTOL, atol=ATOL |
There was a problem hiding this comment.
The variable V_m is not defined in this scope. Based on the context, this should use the existing variable V that is passed to the function, sliced appropriately as V[:, :m].
| k = 1 # Naive arnoldi w/o restart only really works for 1 eigenvalue | ||
| m = m or min(max(2 * k + 1, 20), n) | ||
|
|
||
| V = np.zeros((n, m+1), dtype) |
There was a problem hiding this comment.
Missing parentheses around dtype parameter. Should be V = np.zeros((n, m+1), dtype=dtype) to properly specify the dtype parameter.
| V = np.zeros((n, m+1), dtype) | |
| V = np.zeros((n, m+1), dtype=dtype) |
| V = np.zeros((n, m+1), dtype) | ||
| H = np.zeros((m+1, m), dtype) |
There was a problem hiding this comment.
Missing parentheses around dtype parameter. Should be H = np.zeros((m+1, m), dtype=dtype) to properly specify the dtype parameter.
| V = np.zeros((n, m+1), dtype) | |
| H = np.zeros((m+1, m), dtype) | |
| V = np.zeros((n, m+1), dtype=dtype) | |
| H = np.zeros((m+1, m), dtype=dtype) |
| ritz = RitzDecomposition.from_v_and_h(V, H, k) | ||
| if ritz.approximate_residuals[0] < tol: | ||
| residuals = ritz.compute_true_residuals(A) | ||
| if residuals[0] / max(ritz.values[0], tol) < tol: |
There was a problem hiding this comment.
Division by zero risk if both ritz.values[0] and tol are zero. Consider using max(abs(ritz.values[0]), tol) or adding an explicit check for zero denominator.
| if residuals[0] / max(ritz.values[0], tol) < tol: | |
| if residuals[0] / max(abs(ritz.values[0]), tol) < tol: |
No description provided.