Add LeverageScore distribution - #186
Conversation
Replaces the O(nd^2) single-sketch estimator with the full two-sketch
Algorithm 1 from drineas2012fast: solving M = R^{-1}Pi2 before forming
Omega = A*M avoids materializing A*R^{-1}, giving O(nd log n) instead.
Old single-sketch computation kept as a comment for reference. Also
renames sample_distribution!'s first argument from x to indices to
match the convention used by Uniform, L2Norm, and Agmon (issue #169).
Covers exact mode (Left/Right, including the trivial column-leverage case for tall matrices), approximate mode's argument guards and accuracy, update_distribution!, and sample_distribution!, verified against an independently-derived ground truth rather than the source's own QR call. 100% line coverage on leverage_score.jl.
Wishart bias fix: R'R (from the compressor's sketch) is Wishart-distributed, so naively inverting it overstates every leverage score by an exact factor r1/(r1-d-1). Corrected with the closed-form fix (r1-d-1)/r1, verified numerically to eliminate the bias across several (n, d, r1) regimes. r2 (the second sketch's size) was previously an unprincipled hardcoded formula (20*log(n)) that, worse, was almost always oversized relative to d and thus doing nothing but adding noise. Replaced with an explicit r2 field on LeverageScore/LeverageScoreRecipe, defaulting to nothing (no second sketch, direct A*R^-1); users opt into the two-sketch speedup by setting it, with validation that 1 <= r2 < size(A, 2). Also: renamed the local n_rows/n_cols variables to n/d to stop colliding in meaning with compressor_recipe.n_rows (an unrelated quantity, r1); matched docstring phrasing to the Creates/Updates/Samples convention used by Uniform, L2Norm, and Agmon instead of the inconsistent 'A function that...' phrasing; trimmed several overly long comments and one line over the ~92-char width used elsewhere in the file.
Uniform, L2Norm, and Agmon (written by other collaborators, not renamed by me) all use n_rows/n_cols as their local variable names, never bare n/d. My previous rename to n/d fixed a local reading confusion but broke that established, team-wide convention. Reverted, and resolved the actual point of confusion (compressor_recipe.n_rows meaning something different from A's n_rows) with a one-line comment at r1's definition instead. Also: matched the r2 field's docstring to the same 'description, or nothing in exact mode' pattern used by compressor_recipe, instead of a 'see other struct's docstring' cross-reference with no precedent elsewhere in the module; removed a stray docstring mention of the old r2 = O(log n) formula that no longer applies now that r2 is user-supplied.
r2=nothing previously meant "skip the second sketch, use the full AR^-1 path" unconditionally. It now means "auto-size the second sketch from epsilon and n via the paper's Lemma 1 ε-JLT bound", falling back to the full AR^-1 path only when that bound would not actually be smaller than d. This is what gives approximate mode the paper's O(nd log n) complexity by default instead of requiring users to hand-pick r2 themselves. Explicit r2 values still bypass epsilon entirely. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
drineas2012fast's Lemma 1 is only stated/proven for 0 < epsilon <= 0.5; the previous (0, 1) validation silently accepted epsilon values the JL bound has no guarantee for. Tightened validation and updated the one test that had used epsilon=0.9. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
| factorization of ``A^\\top``. | ||
|
|
||
| If an approximate compressor is provided, leverage scores are estimated using the | ||
| randomized algorithm of [drineas2012fast](@citet): a sketch ``B = \\Pi_1 A`` is |
There was a problem hiding this comment.
For compressors, we have typically used S_1 in stead of \\Pi_1. It might be useful to indicate what the dimensions are (i.e., we should specify r_1.
There was a problem hiding this comment.
S_1 and S_2 instead of Pi; r_1 was defined.
| If an approximate compressor is provided, leverage scores are estimated using the | ||
| randomized algorithm of [drineas2012fast](@citet): a sketch ``B = \\Pi_1 A`` is | ||
| formed via the supplied compressor, and the QR factorization of ``B`` yields ``R``. | ||
| Leverage scores are then the row norms of ``AR^{-1}\\Pi_2``, where |
There was a problem hiding this comment.
We should still use S_2 to indicate a compressor.
| Inverting the sketch-based `R` is a biased estimator of `A'A`'s inverse (`R'R` | ||
| is Wishart-distributed, and matrix inversion is convex, so the naive estimate | ||
| is systematically too large by Jensen's inequality; the exact expectation | ||
| identity is the standard inverse-Wishart mean, e.g. | ||
| [muirhead1982aspects](@citet)). This has an exact, closed-form correction, | ||
| applied internally. It does not fix per-row variance, | ||
| though: rows with small true leverage score are the hardest to pin down to | ||
| tight relative error, since the estimator's noise floor dominates a small | ||
| true value. Approximate mode is best suited to producing a sampling | ||
| distribution (where aggregate weighting matters more than any single row's | ||
| exact value); use exact mode, or a much larger `compressor`, when precision | ||
| matters. |
There was a problem hiding this comment.
I believe this note needs to be clarified a bit. It should indicate what the problem is (R'R derived from a sketch tends to be biased from A'A). If R'R is generated from a Gaussian compressor, then it is known to have a Wishart Distribution, and that there is a bias correction term that can be used. While others compressors may not have the same distribution, the bias correction can still be applied though it may not work.
There was a problem hiding this comment.
Fixed. Note now covers all three: the general problem (R'R biased vs A'A), the Gaussian case (Wishart distribution, exact correction), and the caveat for other S1 types (correction still applied, may not fully resolve the bias).
There was a problem hiding this comment.
There should be three options:
- Compute the leverage scores from QR of A
- Compute the leverage scores from QR of S_1 A, and then using the rows of AR^{-1}
- Compute the leverage scores from QR of S_1 A, and then using the rows of AR^{-1}S_2
There was a problem hiding this comment.
Implemented. compressor2 (S2) is now a field (any Right() compressor, defaults to Gaussian), and LeverageScoreRecipe is parametrized on both compressor-recipe types so each tier gets its own dispatched update_distribution! method. I used a parametric single struct rather than three separately named recipe types (e.g. ExactLSRecipe, ApproxLSRecipe, etc) , to stay consistent with how Agmon handles similar internal complexity (betas), and it extends the same pattern Adrian used in #185 (parameterizing Recipe structs to eliminate dynamic dispatch from abstract field types) instead of introducing a different one. Also SRHTRecipe{C<:Cardinality,...} is already precedent. But I could switch to named subtypes if you'd prefer that instead.
| LeverageScore(; | ||
| cardinality=Undef(), replace=false, compressor=nothing, r2=nothing, | ||
| epsilon=0.5, | ||
| ) |
There was a problem hiding this comment.
This signature should follow the blue style.
Addresses PR #186 review: S1/S2 notation, r1 defined explicitly, clarified Wishart bias-correction note, three tiers now match Vivak's comment, BlueStyle constructor fixed. compressor2 lets any Right() compressor serve as S2 (defaulting to Gaussian), and LeverageScoreRecipe{C1,C2} gives update_distribution! a dispatched, branch-free method per tier instead of runtime nothing-checks. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…cover missing lines in tests
2 participants
Description
Fixes #149
Adds
LeverageScore <: Distributionwith three tiers: exact leverage scores via QR ofA, or approximate ones viaS1alone (AR⁻¹,drineas2012fast), or viaS1and a second compressorS2(AR⁻¹S2) forO(nd log n)cost.S2is a pluggableCompressor(anyRight()-cardinality type), defaulting toGaussiansized by the paper's JL bound. Includes the exact Wishart bias correction for GaussianS1.LeverageScoreRecipeis parametrized on both compressor-recipe types soupdate_distribution!dispatches to a dedicated method per tier.Motivation and Context
Standard sampling primitive for randomized least-squares / CUR / coreset methods.
How has this been tested
285 unit tests (
test/Compressors/Distributions/leverage_score.jl), full suite passes. Docs compiled locally, no errors.Types of changes
Checklists:
Code and Comments
If this PR includes modification to the code base, please select all that apply.
API Documentation
Manual Documentation
Testing
@code_lowered and
@code_typed)