Uh oh!
There was an error while loading. Please reload this page.
Replace OPR's SVDs with Gram-matrix factorizations - #96
Merged
Conversation
Two thirds of an OPR iteration is linear algebra rather than ptychography, and probe_stack grows linearly in frames and modes -- 10.4 GB at 512 px / 890 frames / 6 modes, before the transient the SVD needs on top. On a 32 GB card that is the binding constraint, so this is as much a memory change as a speed one. Both orthogonalization steps are replaced: - orthogonalizeProbeStack ran a full SVD of the (Np^2, nFrames) probe stack, which allocates a (Np^2, nFrames) U it then throws most of away. gram_tsvd gets the same rank-k subspace from the (nFrames, nFrames) Gram matrix instead: 1.8x faster on 5.3x less peak memory at 364/202, 1.9x on 4.3x less at 512/890. - orthogonalizeIncoherentModes ran one small SVD per frame in a Python loop. The batched form factorizes all frames at once: 4.1x at 364/202/4, 1.8x at 512/890/6. The advantage shrinks with size because the loop's cost is mostly per-frame launch overhead, which is what batching removes. It works in frame chunks so the transpose scratch stays bounded rather than allocating a second copy of the whole stack. End to end at 364 px, 202 frames, 4 modes over 30 iterations: 518.8 -> 296.3 ms/iter (1.75x) with peak memory 3.18 -> 1.85 GB (1.72x). Both use svd on the Hermitian Gram matrix rather than eigh. That is mathematically the same -- for a Hermitian positive semi-definite matrix the two coincide, and svd returns descending order for free -- but eigh routes through cupyx.cusolver, which is not importable in every CuPy/CUDA install. On the numerics. The Gram route reaches the same subspace but makes a different arbitrary choice of per-mode global phase than LAPACK, so a naive elementwise comparison of the probe reports relative error 2.0. That is a phase flip, not a different answer. Everything gauge-independent agrees: object 1.4e-07, error metric 4.6e-08, mode powers 4.1e-07 and 1.1e-06, |probe| 8.0e-07, and after aligning the phase the mode vectors agree to 8.2e-07 and 1.7e-04. Over a long run the object does drift, because the solver amplifies any perturbation. At 364/202/30 the Gram route moves the object by 1.6e-02, while perturbing the initial object of the unmodified path by a relative 1e-6 -- below the float32 resolution of the inputs -- moves it by 2.8e-02. The change sits inside the algorithm's own sensitivity to its inputs, and the error metric agrees to 9.3e-05. Defaults change, hence the minor version bump. Set OPR_tsvd_type="numpy" and OPR_fast_orthogonalization=False to recover the previous output; a regression golden pins that path so the escape hatch stays tested. Also fixes orthogonalizeProbeStack to dispatch through getArrayModule instead of calling cp directly, to not mutate the caller's singular values in place on the "numpy" path, and to raise on an unknown OPR_tsvd_type rather than failing later with an unbound name. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Renders the seven derivation steps as LaTeX rather than fenced code blocks, and switches the inline occurrences to match so the note does not carry two notations for the same thing. Also corrects section 7, which cited test names that do not exist: the tests are test_opr_legacy_path_golden and test_opr_gram_matches_legacy_on_gauge_invariant_quantities, in tests/regression/test_opr_regression.py. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for freeto join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Last of the series. Follows #90, #91, #93, #94 and the benchmark PR.
This changes default numerical output — hence
0.2.8→0.3.0rather than a patch bump.Why
Two thirds of an OPR iteration is linear algebra, not ptychography:
orthogonalizeIncoherentModesorthogonalizeProbeStackAnd
probe_stackgrows linearly in frames and modes — 10.4 GB at 512 px / 890 frames /6 modes, before the transient the SVD needs on top. On a 32 GB card that is the binding
constraint, so this is as much a memory change as a speed one.
What changed
orthogonalizeProbeStackran a full SVD of the(Np², nFrames)stack, allocating a(Np², nFrames)U it then throws most of away.gram_tsvdgets the same rank-k subspace fromthe
(nFrames, nFrames)Gram matrix instead.orthogonalizeIncoherentModesran one small SVD per frame in a Python loop. The batchedform factorizes all frames at once, in frame chunks so the transpose scratch stays bounded
rather than allocating a second copy of the whole stack.
orthogonalizeIncoherentModesgram_tsvdvs full SVDThe batched advantage shrinks with size — the loop's cost is mostly per-frame launch overhead,
which is exactly what batching removes; at large sizes the factorization itself dominates.
End to end, 364 px / 202 frames / 4 modes over 30 iterations:
numpyfull SVD (previous default)OPR_tsvd_type="gram"onlyOPR_fast_orthogonalizationonlysvdon the Gram matrix, noteighMathematically identical — for a Hermitian positive semi-definite matrix the SVD and the
eigendecomposition coincide, and
svdreturns descending order for free. The practical reason isthat
cupy.linalg.eighis not importable in every CuPy/CUDA install (it routes throughcupyx.cusolver; here it fails onlibcusolver.so.11). The implementation deliberately does notdepend on it.
On the numerics — please read this bit
A naive elementwise comparison of the probe reports relative error 2.0, which looks
catastrophic and is not. Mode vectors are defined only up to a per-mode global phase, and the
Gram route makes a different arbitrary choice than LAPACK. A phase flip is exactly relative
error 2.
Everything gauge-independent agrees:
|probe|Over a longer run the object does drift, because the solver amplifies any perturbation. At
364 px / 202 frames / 30 iterations the Gram route moves the object by 1.6e-02 — while
perturbing the initial object of the unmodified path by a relative 1e-6, below the float32
resolution of the inputs, moves it by 2.8e-02. The change sits inside the algorithm's own
sensitivity to its inputs, and the error metric agrees to 9.3e-05.
Escape hatch, and it is tested
recovers the previous output. A regression golden pins that path, so the claim is verified rather
than asserted.
Tests
tests/Engines/test_opr_linalg.py—gram_tsvdunit tests, on CPU so CI covers them eventhough OPR itself is GPU-only. Includes rank clamping and rank-deficient input.
tests/regression/test_opr_regression.py— the legacy golden, and the gauge-invariantcomparison above.
Suite goes 61 → 71 passing.
Incidental fixes in the same functions
orthogonalizeProbeStackdispatches throughgetArrayModuleinstead of callingcpdirectly"numpy"path no longer mutates the caller's singular values in placeOPR_tsvd_typenow raises, instead of failing later with an unbound nameNote
notes/gram_orthogonalization.mdis a full walkthrough of why the Gram route gives the sameanswer and where the speedup comes from — written for someone who has not thought about this
before.