From 8706f0674eb87e01067ad0f43b3f201b41e4cad5 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 24 Aug 2024 23:12:27 -0400 Subject: [PATCH] a lot --- RandBLAS/base.hh | 53 ++++++-- RandBLAS/dense_skops.hh | 65 ++++------ RandBLAS/sparse_skops.hh | 69 ++++++----- .../total-least-squares/tls_sparse_skop.cc | 12 +- rtd/source/api_reference/skops_and_dists.rst | 104 ++++++++-------- rtd/source/tutorial/index.rst | 2 +- rtd/source/tutorial/sketch_updates.rst | 117 ++++++++++++++++++ rtd/source/updates/index.rst | 2 +- .../randblas_rtd/static/theme_overrides.css | 58 +++------ test/test_datastructures/test_sparseskop.cc | 10 +- .../test_spmats/test_coo.cc | 5 +- test/test_handrolled_lapack.cc | 6 +- test/test_matmul_cores/test_lskges.cc | 29 ++--- test/test_matmul_cores/test_rskges.cc | 15 +-- 14 files changed, 316 insertions(+), 231 deletions(-) create mode 100644 rtd/source/tutorial/sketch_updates.rst diff --git a/RandBLAS/base.hh b/RandBLAS/base.hh index 4c2cf351..0cf9138b 100644 --- a/RandBLAS/base.hh +++ b/RandBLAS/base.hh @@ -283,10 +283,33 @@ enum class MajorAxis : char { // ============================================================================= /// @verbatim embed:rst:leading-slashes /// -/// TODO: take the expectation properties in scaled partial isometries docs and move here. -/// (Or, rather, move above.) +/// **Mathematical description** /// -/// Any object :math:`\ttt{D}` of type :math:`\ttt{SkDist}` has the following attributes. +/// Matrices sampled from sketching distributions in RandBLAS are mean-zero +/// and have covariance matrices that are proportional to the identity. That is, +/// if :math:`\D` is a distribution over :math:`r \times c` matrices and +/// :math:`\mtxS` is a sample from :math:`\D`, then +/// :math:`\mathbb{E}\mtxS = \mathbf{0}_{r \times c}` and +/// +/// .. math:: +/// :nowrap: +/// +/// \begin{gather} +/// \alpha^2 \cdot \mathbb{E}\left[ \mtxS^T\mtxS \right]=\mathbf{I}_{c \times c}& \nonumber \\ +/// \,\beta^2 \cdot \mathbb{E}\left[ \mtxS{\mtxS}^T\, \right]=\mathbf{I}_{r \times r}& \nonumber +/// \end{gather} +/// +/// hold for some :math:`\alpha > 0` and :math:`\beta > 0`. +/// +/// The *isometry scale* of the distribution +/// is :math:`\gamma := \alpha` if :math:`c \geq r` and :math:`\gamma := \beta` otherwise. If you want to +/// sketch in a way that preserves squared norms in expectation, then you should sketch with +/// a scaled sample :math:`\gamma \mtxS` rather than the sample itself. +/// +/// **Programmatic description** +/// +/// A variable :math:`\ttt{D}` of a type that conforms to the +/// :math:`\ttt{SketchingDistribution}` concept has the following attributes. /// /// .. list-table:: /// :widths: 25 30 40 @@ -303,11 +326,13 @@ enum class MajorAxis : char { /// - samples from :math:`\ttt{D}` have this many columns /// * - :math:`\ttt{D.major_axis}` /// - :math:`\ttt{const MajorAxis}` -/// - Indicate the nature of statistical independence among a sample's entries. +/// - Implementation-dependent; see MajorAxis documentation. +/// * - :math:`\ttt{D.isometry_scale}` +/// - :math:`\ttt{const double}` +/// - See above. /// -/// -/// Note that we only allow one major axis. Although it's not obvious, enforcing a single major axis -/// is relevant even for distributions over matrices with i.i.d. entries. +/// Note that the isometry scale is always stored in double precision; this has no bearing +/// on the precision of sketching operators that are sampled from a :math:`\ttt{SketchingDistribution}`. /// /// @endverbatim template @@ -315,7 +340,14 @@ concept SketchingDistribution = requires(SkDist D) { { D.n_rows } -> std::same_as; { D.n_cols } -> std::same_as; { D.major_axis } -> std::same_as; + { D.isometry_scale } -> std::same_as; }; +#else +/// .. math:: +/// +/// {\D}\texttt{.isometry\_scale} = \begin{cases} \alpha &\text{ if } r \leq c \\ \beta &\text{ if } r > c \end{cases}~. +#define SketchingDistribution typename +#endif // ============================================================================= /// \fn isometry_scale_factor(SkDist D) @@ -348,9 +380,7 @@ concept SketchingDistribution = requires(SkDist D) { /// @endverbatim template inline T isometry_scale_factor(SkDist D); -#else -#define SketchingDistribution typename -#endif + #ifdef __cpp_concepts // ============================================================================= @@ -365,8 +395,9 @@ template concept SketchingOperator = requires(SKOP S) { { S.n_rows } -> std::same_as; { S.n_cols } -> std::same_as; + { S.dist } -> SketchingDistribution; { S.seed_state } -> std::same_as; - { S.seed_state } -> std::same_as; + { S.next_state } -> std::same_as; }; #else #define SketchingOperator typename diff --git a/RandBLAS/dense_skops.hh b/RandBLAS/dense_skops.hh index 52ac7b8c..0740d658 100644 --- a/RandBLAS/dense_skops.hh +++ b/RandBLAS/dense_skops.hh @@ -190,6 +190,11 @@ RNGState compute_next_state(DD dist, RNGState state) { return state; } +template +inline double isometry_scale(DDN dn, int64_t n_rows, int64_t n_cols) { + return (dn == DDN::BlackBox) ? 1.0 : std::pow(std::min(n_rows, n_cols), -0.5); +} + } // end namespace RandBLAS::dense @@ -228,45 +233,25 @@ struct DenseDist { /// Matrices drawn from this distribution have this many columns. const int64_t n_cols; - // --------------------------------------------------------------------------- - /// The distribution used for the entries of the sketching operator. - const DenseDistName family; - // --------------------------------------------------------------------------- /// This member affects whether samples from this distribution have their /// entries filled row-wise or column-wise. While there is no statistical /// difference between these two filling orders, there are situations /// where one order or the other might be preferred. - /// /// @verbatim embed:rst:leading-slashes - /// .. dropdown:: *Notes for experts* - /// :animate: fade-in-slide-down - /// - /// Deciding the value of this member is only needed - /// in algorithms where (1) there's a need to iteratively generate panels of - /// a larger sketching operator and (2) one of larger operator's dimensions - /// cannot be known before the iterative process starts. - /// - /// Essentially, a column-wise fill order lets us stack operators horizontally - /// in a consistent way, while row-wise fill order lets us stack vertically - /// in a consistent way. The mapping from major_axis to fill order is given below. - /// - /// .. list-table:: - /// :widths: 34 33 33 - /// :header-rows: 1 - /// - /// * - - /// - :math:`\texttt{major_axis} = \texttt{Long}` - /// - :math:`\texttt{major_axis} = \texttt{Short}` - /// * - :math:`\texttt{n_rows} > \texttt{n_cols}` - /// - column-wise - /// - row-wise - /// * - :math:`\texttt{n_rows} \leq \texttt{n_cols}` - /// - row-wise - /// - column-wise + /// See our tutorial section on :ref:`updating sketches ` for more information. /// @endverbatim const MajorAxis major_axis; + // --------------------------------------------------------------------------- + /// A sketching operator sampled from this distribution should be multiplied + /// by this constant in order for sketching to preserve norms in expectation. + const double isometry_scale; + + // --------------------------------------------------------------------------- + /// The distribution used for the entries of the sketching operator. + const DenseDistName family; + // --------------------------------------------------------------------------- /// A distribution over matrices of shape (n_rows, n_cols) with entries drawn /// iid from either the default choice of standard normal distribution, or from @@ -276,14 +261,17 @@ struct DenseDist { int64_t n_rows, int64_t n_cols, DenseDistName dn = DenseDistName::Gaussian - ) : n_rows(n_rows), n_cols(n_cols), family(dn), major_axis( (dn == DenseDistName::BlackBox) ? MajorAxis::Undefined : MajorAxis::Long) { } + ) : n_rows(n_rows), n_cols(n_cols), + major_axis( (dn == DenseDistName::BlackBox) ? MajorAxis::Undefined : MajorAxis::Long), + isometry_scale(dense::isometry_scale(dn, n_rows, n_cols)), family(dn) { } DenseDist( int64_t n_rows, int64_t n_cols, DenseDistName dn, MajorAxis ma - ) : n_rows(n_rows), n_cols(n_cols), family(dn), major_axis(ma) { + ) : n_rows(n_rows), n_cols(n_cols), major_axis(ma), + isometry_scale(dense::isometry_scale(dn, n_rows, n_cols)), family(dn) { if (dn == DenseDistName::BlackBox) { randblas_require(ma == MajorAxis::Undefined); } else { @@ -315,16 +303,6 @@ inline int64_t major_axis_length(const DenseDist &D) { std::max(D.n_rows, D.n_cols) : std::min(D.n_rows, D.n_cols); } -template -inline T isometry_scale_factor(DenseDist D) { - if (D.family == DenseDistName::BlackBox) { - throw std::runtime_error("Unrecognized distribution."); - } - // When we sample from the scalar distributions we always - // scale things so they're variance-1. - return std::pow((T) std::min(D.n_rows, D.n_cols), -0.5); -} - // ============================================================================= /// A sample from a distribution over dense sketching operators. @@ -341,7 +319,7 @@ struct DenseSkOp { // --------------------------------------------------------------------------- /// The distribution from which this sketching operator is sampled. /// This member specifies the number of rows and columns of the sketching - /// operator. + /// operator, as well as whether it is generated row-wise or column-wise. const DenseDist dist; // --------------------------------------------------------------------------- @@ -418,6 +396,7 @@ struct DenseSkOp { } }; +static_assert(SketchingDistribution); // ============================================================================= /// @verbatim embed:rst:leading-slashes diff --git a/RandBLAS/sparse_skops.hh b/RandBLAS/sparse_skops.hh index e3cba94b..4c356999 100644 --- a/RandBLAS/sparse_skops.hh +++ b/RandBLAS/sparse_skops.hh @@ -105,6 +105,18 @@ static RNGState repeated_fisher_yates( return RNGState {ctr, key}; } +inline double isometry_scale(MajorAxis ma, int64_t vec_nnz, int64_t n_rows, int64_t n_cols) { + if (ma == MajorAxis::Short) { + return std::pow(vec_nnz, -0.5); + } else if (ma == MajorAxis::Long) { + double minor_ax_len = std::min(n_rows, n_cols); + double major_ax_len = std::max(n_rows, n_cols); + return std::sqrt( major_ax_len / (vec_nnz * minor_ax_len) ); + } else { + throw std::invalid_argument("Cannot compute the isometry scale for a sparse operator with unspecified major axis."); + } +} + } namespace RandBLAS { @@ -121,6 +133,20 @@ struct SparseDist { /// Matrices drawn from this distribution have this many columns. const int64_t n_cols; + // --------------------------------------------------------------------------- + /// Constrains the sparsity pattern of matrices drawn from this distribution. + /// + /// Having major_axis==Short results in sketches are more likely to contain + /// useful geometric information, without making assumptions about the data + /// being sketched. + /// + const MajorAxis major_axis = MajorAxis::Short; + + // --------------------------------------------------------------------------- + /// A sketching operator sampled from this distribution should be multiplied + /// by this constant in order for sketching to preserve norms in expectation. + const double isometry_scale; + // --------------------------------------------------------------------------- /// If this distribution is short-axis major, then matrices sampled from /// it will have exactly \math{\texttt{vec_nnz}} nonzeros per short-axis @@ -134,14 +160,15 @@ struct SparseDist { /// const int64_t vec_nnz; - // --------------------------------------------------------------------------- - /// Constrains the sparsity pattern of matrices drawn from this distribution. - /// - /// Having major_axis==Short results in sketches are more likely to contain - /// useful geometric information, without making assumptions about the data - /// being sketched. - /// - const MajorAxis major_axis = MajorAxis::Short; + SparseDist( + int64_t n_rows, + int64_t n_cols, + MajorAxis ma, + int64_t vec_nnz + ) : n_rows(n_rows), n_cols(n_cols), major_axis(ma), + isometry_scale(sparse::isometry_scale(ma, vec_nnz, n_rows, n_cols)), + vec_nnz(vec_nnz) {} + }; template @@ -151,18 +178,6 @@ inline RNGState repeated_fisher_yates( return sparse::repeated_fisher_yates(state, k, n, r, indices, (sint_t*) nullptr, (double*) nullptr); } -template -inline T isometry_scale_factor(SparseDist D) { - T vec_nnz = (T) D.vec_nnz; - if (D.major_axis == MajorAxis::Short) { - return std::pow(vec_nnz, -0.5); - } else { - T minor_ax_len = (T) std::min(D.n_rows, D.n_cols); - T major_ax_len = (T) std::max(D.n_rows, D.n_cols); - return std::sqrt( major_ax_len / (vec_nnz * minor_ax_len) ); - } -} - template RNGState compute_next_state(SparseDist dist, RNGState state) { int64_t minor_len; @@ -176,7 +191,6 @@ RNGState compute_next_state(SparseDist dist, RNGState state) { return state; } - // ============================================================================= /// A sample from a prescribed distribution over sparse matrices. /// @@ -193,7 +207,7 @@ struct SparseSkOp { // --------------------------------------------------------------------------- /// The distribution from which this sketching operator is sampled. /// This member specifies the number of rows and columns of the sketching - /// operator. + /// operator. It also controls the sparsity pattern and the sparsity level. const SparseDist dist; // --------------------------------------------------------------------------- @@ -376,6 +390,8 @@ struct SparseSkOp { } }; +static_assert(SketchingDistribution); + // ============================================================================= /// Performs the work in sampling S from its underlying distribution. This /// entails populating S.rows, S.cols, and S.vals with COO-format sparse matrix @@ -501,12 +517,9 @@ COOMatrix coo_view_of_skop(SparseSkOp &S) { template static auto transpose(SparseSkOp const &S) { randblas_require(S.known_filled); - SparseDist dist = { - .n_rows = S.dist.n_cols, - .n_cols = S.dist.n_rows, - .vec_nnz = S.dist.vec_nnz, - .major_axis = S.dist.major_axis - }; + SparseDist dist( + S.dist.n_cols, S.dist.n_rows, S.dist.major_axis, S.dist.vec_nnz + ); SparseSkOp St(dist, S.seed_state, S.cols, S.rows, S.vals); St.next_state = S.next_state; return St; diff --git a/examples/total-least-squares/tls_sparse_skop.cc b/examples/total-least-squares/tls_sparse_skop.cc index 44c5c3e6..bca45e05 100644 --- a/examples/total-least-squares/tls_sparse_skop.cc +++ b/examples/total-least-squares/tls_sparse_skop.cc @@ -139,12 +139,12 @@ int main(int argc, char* argv[]){ // Sample the sketching operator auto time_constructsketch1 = high_resolution_clock::now(); - RandBLAS::SparseDist Dist = { - .n_rows = sk_dim, // Number of rows of the sketching operator - .n_cols = m, // Number of columns of the sketching operator - .vec_nnz = 8, // Number of non-zero entires per major-axis vector - .major_axis = RandBLAS::MajorAxis::Short // A "SASO" (aka SJLT, aka OSNAP, aka generalized CountSketch) - }; + RandBLAS::SparseDist Dist( + sk_dim, // Number of rows of the sketching operator + m, // Number of columns of the sketching operator + RandBLAS::MajorAxis::Short, // A "SASO" (aka SJLT, aka OSNAP, aka generalized CountSketch) + 8 // Number of non-zero entires per column + ); uint32_t seed = 1997; RandBLAS::SparseSkOp S(Dist, seed); RandBLAS::fill_sparse(S); diff --git a/rtd/source/api_reference/skops_and_dists.rst b/rtd/source/api_reference/skops_and_dists.rst index 0afdd53c..2393f17f 100644 --- a/rtd/source/api_reference/skops_and_dists.rst +++ b/rtd/source/api_reference/skops_and_dists.rst @@ -17,62 +17,55 @@ .. |mtxy| mathmacro:: \mathbf{y} .. |ttt| mathmacro:: \texttt -****************************************************** -Fundamentals of sketching and random number generation -****************************************************** +******************************************************************** +Fundamentals +******************************************************************** -.. _rngstate_api: + .. very similar effects can be acheived in C, Fortran, or Julia. While there are + .. certainly some popular programming languages that don't support this kind of API + .. (e.g., MATLAB, Python, and R), accessing RandBLAS from these languages should + .. be mediated operator-overloaded objects in a way that's analogous to how one + .. would access BLAS. +RandBLAS has a polymorphic free-function API. We have spent a significant amount of +effort on minimizing the number of RandBLAS-specific datastructures needed in order +to acheive the polymorphic API. -Sketching in RandBLAS concerns linear maps (sketching operators) that take vectors in high-dimensional -coordinate representations to vectors in low-dimensional coordinate representations. +RandBLAS is very light on C++ idioms. The main C++ idioms we use are +templating and function overloading, plus some mild memory management +with destructor methods for structs. -We assume that the initial vectors are given as the rows or columns of some matrix. -This lets us describe sketching operators as matrices that are wide (having more columns than rows) -or tall (having more rows than columns). -Square sketching operators are not forbidden, but they are certainly not our focus. +:ref:`Test of sketch updates `. + + +Abstractions +============ .. dropdown:: Distributions over random matrices :animate: fade-in-slide-down :color: light + :open: .. doxygenconcept:: RandBLAS::SketchingDistribution :project: RandBLAS - .. dropdown:: Scaling to obtain partial isometries - :animate: fade-in-slide-down - :color: light - - .. doxygenfunction:: RandBLAS::isometry_scale_factor(SkDist D) - :project: RandBLAS - - .. dropdown:: Details on MajorAxis - :animate: fade-in-slide-down - :color: light - - .. doxygenenum:: RandBLAS::MajorAxis - :project: RandBLAS - -.. dropdown:: Sketching operators - :animate: fade-in-slide-down - :color: light - - .. doxygenconcept:: RandBLAS::SketchingOperator +.. dropdown:: Details on MajorAxis + :animate: fade-in-slide-down + :color: light + + .. doxygenenum:: RandBLAS::MajorAxis :project: RandBLAS -.. dropdown:: States of random number generators +.. dropdown:: Shared aspects of the sketching operator interface :animate: fade-in-slide-down :color: light - - .. doxygenstruct:: RandBLAS::RNGState + + .. doxygenconcept:: RandBLAS::SketchingOperator :project: RandBLAS - :members: - -.. _densedist_and_denseskop_api: -Dense sketching: Gaussians et al. -================================= +Distributions +============= .. dropdown:: DenseDist : a distribution over matrices with i.i.d., mean-zero, variance-one entries :animate: fade-in-slide-down @@ -85,34 +78,39 @@ Dense sketching: Gaussians et al. .. doxygenenum:: RandBLAS::DenseDistName :project: RandBLAS -.. dropdown:: DenseSkOp : a sample from a DenseDist +.. dropdown:: SparseDist : a distribution over structured sparse matrices :animate: fade-in-slide-down :color: light - .. doxygenstruct:: RandBLAS::DenseSkOp + .. doxygenstruct:: RandBLAS::SparseDist :project: RandBLAS - :members: + :members: - *Memory management* - .. doxygenfunction:: RandBLAS::fill_dense(DenseSkOp &S) - :project: RandBLAS +Random states and sketching operators +===================================== - .. doxygenfunction:: RandBLAS::fill_dense(blas::Layout layout, const DenseDist &D, int64_t n_rows, int64_t n_cols, int64_t S_ro, int64_t S_co, T *buff, const RNGState &seed) - :project: RandBLAS - -.. _sparsedist_and_sparseskop_api: +.. dropdown:: RNGState + :animate: fade-in-slide-down + :color: light -Sparse sketching: CountSketch et al. -==================================== + .. doxygenstruct:: RandBLAS::RNGState + :project: RandBLAS + :members: -.. dropdown:: SparseDist : a distribution over structured sparse matrices +.. dropdown:: DenseSkOp : a sample from a DenseDist :animate: fade-in-slide-down :color: light - .. doxygenstruct:: RandBLAS::SparseDist + .. doxygenstruct:: RandBLAS::DenseSkOp :project: RandBLAS - :members: + :members: + + .. doxygenfunction:: RandBLAS::fill_dense(DenseSkOp &S) + :project: RandBLAS + + .. doxygenfunction:: RandBLAS::fill_dense(blas::Layout layout, const DenseDist &D, int64_t n_rows, int64_t n_cols, int64_t S_ro, int64_t S_co, T *buff, const RNGState &seed) + :project: RandBLAS .. dropdown:: SparseSkOp : a sample from a SparseDist :animate: fade-in-slide-down @@ -124,5 +122,3 @@ Sparse sketching: CountSketch et al. .. doxygenfunction:: RandBLAS::fill_sparse(SparseSkOp &S) :project: RandBLAS - - diff --git a/rtd/source/tutorial/index.rst b/rtd/source/tutorial/index.rst index aa0662fb..e1aa0677 100644 --- a/rtd/source/tutorial/index.rst +++ b/rtd/source/tutorial/index.rst @@ -51,6 +51,6 @@ It even lets you compute products against *submatrices* of sketching operators w Background on GEMM Defining a sketching distribution Sampling a sketching operator - Updating a sketch + Updating a sketch The meaning of "submat(・)" in RandBLAS documentation diff --git a/rtd/source/tutorial/sketch_updates.rst b/rtd/source/tutorial/sketch_updates.rst new file mode 100644 index 00000000..d6f89ba1 --- /dev/null +++ b/rtd/source/tutorial/sketch_updates.rst @@ -0,0 +1,117 @@ + + + .. |denseskop| mathmacro:: \texttt{DenseSkOp} + .. |seedstate| mathmacro:: \texttt{seed_state} + .. |nextstate| mathmacro:: \texttt{next_state} + .. |mtx| mathmacro:: {} + +.. _sketch_updates: + +********************************************************************************************* +Updating sketches with dense sketching operators +********************************************************************************************* + +RandBLAS makes it easy to *implicitly* extend an initial sketching +operator :math:`\mtx{S}_1` into an augmented operator :math:`\mtx{S} = [\mtx{S}_1; \mtx{S}_2]` or :math:`\mtx{S} = [\mtx{S}_1, \mtx{S}_2]`. +There are four scenarios that you can find yourself in where +this can be done without generating :math:`S` from scratch. +In all four scenarios, the idea is to +use the :math:`\nextstate` of :math:`\mtx{S}_1` as the +:math:`\seedstate` of :math:`\mtx{S}_2`. + +There are two reasons why you'd want to +extend a sketching operator; you might be trying to improve statistical +quality by increasing sketch size, or you might be +incorporating new data into a sketch of fixed size. +The unifying perspective on the former scenarios is that they both add +*long-axis vectors* to the sketching operator. +The unifying perspective on +Scenarios 3 and 4 is that they both add *short-axis vectors* to the +sketching operator. + +:math:`\texttt{DenseDist}` objects have a :math:`\texttt{major_axis}` member, which states +whether operators sampled from that distribution are short-axis or +long-axis major. So when you specify the major axis for a sketching +operator, you're basically saying whether you want to keep open the possibility of +improving the statistical quality of a sketch or updating a sketch to +incorporate more data. + + +Increase sketch size by adding long-axis vectors. +================================================= + + Suppose :math:`\mtx{S}_1` is a *wide* :math:`d_1 \times m` row-wise + :math:`\denseskop` with seed :math:`c`. It's easy to generate a + :math:`d_2\times m` row-wise :math:`\denseskop` :math:`\mtx{S}_2` in such a way that + :math:`\mtx{S} = [\mtx{S}_1; \mtx{S}_2]` is the same as the :math:`(d_1 + d_2) \times m` row-wise + :math:`\denseskop` with seed :math:`c`. + + This scenario arises if we have a fixed data matrix :math:`\mtx{A}`, an initial + sketch :math:`\mtx{B}_1 = \mtx{S}_1 \mtx{A}`, and we decide we want a larger sketch for + statistical reasons. The updated sketch :math:`\mtx{B} = \mtx{S} \mtx{A}` can be expressed as + + .. math:: + + \mtx{B} = \begin{bmatrix} \mtx{S}_1 \mtx{A} \\ \mtx{S}_2 \mtx{A} \end{bmatrix}. + + Now suppose :math:`\mtx{S}_1` a *tall* :math:`n \times d_1` column-wise :math:`\denseskop` + with seed :math:`c`. We can easily generate an :math:`n\times d_2` column-wise + :math:`\denseskop` :math:`\mtx{S}_2` so that :math:`\mtx{S} = [\mtx{S}_1, \mtx{S}_2]` is the same + as the :math:`d \times (n_1 + n_2)` column-wise :math:`\denseskop` with seed :math:`c`. + + This arises we have a fixed data matrix :math:`\mtx{A}`, an initial sketch :math:`\mtx{B}_1 = \mtx{A} \mtx{S}_1`, + and we decide we want a larger sketch for statistical reasons. The + updated sketch :math:`\mtx{B} = \mtx{A}\mtx{S}` can be expressed as + + .. math:: + + \mtx{B} = \begin{bmatrix} \mtx{A} \mtx{S}_1 & \mtx{A} \mtx{S}_2 \end{bmatrix}. + +If :math:`(\mtx{S}_1, \mtx{S}_2, \mtx{S})` satisfy the assumptions above, then the final sketch +:math:`B` will be the same regardless of whether we computed the sketch in one +step or two steps. This is desirable for benchmarking and debugging +RandNLA algorithms where the sketch size is a tuning parameter. + + +Accomodate more data by adding short-axis vectors. +================================================== + + Suppose :math:`\mtx{S}_1` is a *tall* :math:`n_1 \times d` row-wise + :math:`\denseskop` with seed :math:`c`. It's easy to generate an :math:`n_2\times d` + row-wise :math:`\denseskop` :math:`\mtx{S}_2` in such a way that + :math:`\mtx{S} = [\mtx{S}_1; \mtx{S}_2]` is the same as the :math:`(n_1 + n_2) \times d` row-wise + :math:`\denseskop` with seed :math:`c`. + + This situation arises if we have an initial data matrix :math:`\mtx{A}_1`, an initial sketch + :math:`\mtx{B}_1 = \mtx{A}_1 \mtx{S}_1`, and then a new matrix :math:`\mtx{A}_2` arrives in such a way that we + want a sketch of :math:`\mtx{A} = [\mtx{A}_1, \mtx{A}_2]`. To compute :math:`\mtx{B} = \mtx{A}\mtx{S}`, we update :math:`\mtx{B}_1` + according to the formula + + .. math:: + + \mtx{B} = \begin{bmatrix} \mtx{A}_1 & \mtx{A}_2 \end{bmatrix} \begin{bmatrix} \mtx{S}_1 \\ \mtx{S}_2 \end{bmatrix} = \mtx{B}_1 + \mtx{A}_2 \mtx{S}_2. + + Now, suppose instead :math:`\mtx{S}_1` is a *wide* :math:`d \times m_1` column-wise + :math:`\denseskop` with seed :math:`c`. It's easy to generate a + :math:`d \times m_2` column-wise :math:`\denseskop` :math:`\mtx{S}_2` so that + :math:`\mtx{S} = [\mtx{S}_1, \mtx{S}_2]` is the same as the :math:`d \times (m_1 + m_2)` column-wise + :math:`\denseskop` with seed :math:`c`. + + This situation arises if we have an initial data matrix :math:`\mtx{A}_1`, an + initial sketch :math:`\mtx{B}_1 = \mtx{S}_1 \mtx{A}_1`, and then a new matrix :math:`\mtx{A}_2` arrives in + such a way that we want a sketch of :math:`A = [\mtx{A}_1; \mtx{A}_2]`. To compute :math:`\mtx{B} = \mtx{S}\mtx{A}`, + we update :math:`\mtx{B}_1` according to the formula + + .. math:: + + \mtx{B} = \begin{bmatrix} \mtx{S}_1 & \mtx{S}_2 \end{bmatrix} \begin{bmatrix} \mtx{A}_1 \\ \mtx{A}_2 \end{bmatrix} = \mtx{B}_1 + \mtx{S}_2 \mtx{A}_2. + +If :math:`(\mtx{S}_1, \mtx{S}_2, \mtx{S})` satisfy the assumptions above, then :math:`\mtx{B}` will be the +same as though we started with all of :math:`\mtx{A}` from the very beginning. This +is useful for benchmarking and debugging RandNLA algorithms that involve +operating on data matrices that increase in size over some number of iterations. + + +Porque no los dos? Work with giant, implicit operators. +========================================================== + diff --git a/rtd/source/updates/index.rst b/rtd/source/updates/index.rst index d8a4667a..236685ea 100644 --- a/rtd/source/updates/index.rst +++ b/rtd/source/updates/index.rst @@ -8,7 +8,7 @@ RandBLAS upon request, no matter how old. With any luck, RandBLAS will grow enou in the future that we will change this policy to support a handful of versions at a time. -RandBLAS follows `Semantic Versioning `_. +RandBLAS follows `Semantic Versioning `_. RandBLAS 0.2 diff --git a/rtd/themes/randblas_rtd/static/theme_overrides.css b/rtd/themes/randblas_rtd/static/theme_overrides.css index 3e78816d..c6141883 100644 --- a/rtd/themes/randblas_rtd/static/theme_overrides.css +++ b/rtd/themes/randblas_rtd/static/theme_overrides.css @@ -1,18 +1,14 @@ @import 'css/theme.css'; /* override table width restrictions */ -/* .wy-table-responsive table td, .wy-table-responsive table th { +.wy-table-responsive table td, .wy-table-responsive table th { white-space: normal; -} */ - -.breatheparameterlist li tt + p { - display: inline; - } +} .wy-table-responsive { margin-bottom: 24px; max-width: 100%; - overflow: scroll; + overflow: auto; } .wy-nav-side { @@ -37,8 +33,8 @@ padding: 1.618em 3.236em; height: 100%; max-width: 1000px; - margin: auto; - /* text-align: justify; */ + margin-left: 0; + text-align: justify; } .wy-nav-content-wrap.collapsed { @@ -49,52 +45,28 @@ .wy-nav-content.collapsed { margin-left: 0; - max-width: none; + max-width: 1000px; margin: 0; width: 100%; } .math { text-align: left; - overflow: scroll; + overflow: auto; } .eqno { float: right; } - -/* .cpp.class, .cpp.class dl, .cpp.class dd, .cpp.class dt { - margin-left: 0.5em; - padding-left: 0; -} - -.cpp.enumeration, .cpp.enumeration dl, .cpp.enumeration dd, .cpp.enumeration dt { - margin-left: 0.5em; - padding-left: 0; -} */ - +/* +The standard rendering of c++ doyxgen comments involves large margins for many elements. +The styles below set override the margin choice for all (?) relevant HTML elements. +*/ +.cpp.function, .cpp.function dl, .cpp.function dd, .cpp.function dt, +.cpp.struct, .cpp.struct dl, .cpp.struct dd, .cpp.struct dt, +.cpp.concept, .cpp.concept dl, .cpp.concept dd, .cpp.concept dt, +.cpp.var, .cpp.var dl, .cpp.var dd, .cpp.var dt, .cpp.enum-class, .cpp.enum-class dl, .cpp.enum-class dd, .cpp.enum-class dt { margin-left: 0.5em; - /* padding-left: 0; */ -} - -.cpp.var, .cpp.var dl, .cpp.var dd, .cpp.var dt { - margin-left: 0.5em; - /* padding-left: 0; */ -} - -.cpp.concept, .cpp.concept dl, .cpp.concept dd, .cpp.concept dt { - margin-left: 0.5em; - /* padding-left: 0; */ -} - -.cpp.struct, .cpp.struct dl, .cpp.struct dd, .cpp.struct dt { - margin-left: 0.5em; - /* padding-left: 0; */ -} - -.cpp.function, .cpp.function dl, .cpp.function dd, .cpp.function dt { - margin-left: 0.5em; - /* padding-left: 0; */ } diff --git a/test/test_datastructures/test_sparseskop.cc b/test/test_datastructures/test_sparseskop.cc index 7709de52..cb31a5b5 100644 --- a/test/test_datastructures/test_sparseskop.cc +++ b/test/test_datastructures/test_sparseskop.cc @@ -78,9 +78,8 @@ class TestSparseSkOpConstruction : public ::testing::Test template void proper_saso_construction(int64_t d, int64_t m, int64_t key_index, int64_t nnz_index) { using RNG = RandBLAS::SparseSkOp::state_t::generator; - RandBLAS::SparseSkOp S0( - {d, m, vec_nnzs[nnz_index], RandBLAS::MajorAxis::Short}, keys[key_index] - ); + RandBLAS::SparseDist D0(d, m, RandBLAS::MajorAxis::Short, vec_nnzs[nnz_index]); + RandBLAS::SparseSkOp S0(D0, keys[key_index]); RandBLAS::fill_sparse(S0); if (d < m) { check_fixed_nnz_per_col(S0); @@ -92,9 +91,8 @@ class TestSparseSkOpConstruction : public ::testing::Test template void proper_laso_construction(int64_t d, int64_t m, int64_t key_index, int64_t nnz_index) { using RNG = RandBLAS::SparseSkOp::state_t::generator; - RandBLAS::SparseSkOp S0( - {d, m, vec_nnzs[nnz_index], RandBLAS::MajorAxis::Long}, keys[key_index] - ); + RandBLAS::SparseDist D0(d, m, RandBLAS::MajorAxis::Long, vec_nnzs[nnz_index]); + RandBLAS::SparseSkOp S0(D0, keys[key_index]); RandBLAS::fill_sparse(S0); if (d < m) { check_fixed_nnz_per_row(S0); diff --git a/test/test_datastructures/test_spmats/test_coo.cc b/test/test_datastructures/test_spmats/test_coo.cc index a9bb0bbd..03ca0e28 100644 --- a/test/test_datastructures/test_spmats/test_coo.cc +++ b/test/test_datastructures/test_spmats/test_coo.cc @@ -139,9 +139,8 @@ class Test_SkOp_to_COO : public ::testing::Test { template void sparse_skop_to_coo(int64_t d, int64_t m, int64_t key_index, int64_t nnz_index, RandBLAS::MajorAxis ma) { - RandBLAS::SparseSkOp S( - {d, m, vec_nnzs[nnz_index], ma}, keys[key_index] - ); + RandBLAS::SparseDist D(d, m, ma, vec_nnzs[nnz_index]); + RandBLAS::SparseSkOp S(D, keys[key_index]); auto A = RandBLAS::sparse::coo_view_of_skop(S); EXPECT_EQ(S.dist.n_rows, A.n_rows); diff --git a/test/test_handrolled_lapack.cc b/test/test_handrolled_lapack.cc index 21f4fa1f..3a97d767 100644 --- a/test/test_handrolled_lapack.cc +++ b/test/test_handrolled_lapack.cc @@ -32,7 +32,7 @@ class TestHandrolledCholesky : public ::testing::Test { DenseDist D(m, n); std::vector A(n*n); std::vector B(m*n); - T iso_scale = std::pow(RandBLAS::isometry_scale_factor(D), 2); + T iso_scale = std::pow(D.isometry_scale, 2); RNGState state(key); RandBLAS::fill_dense(D, B.data(), state); std::vector C(B); @@ -177,7 +177,7 @@ class TestHandrolledQR : public ::testing::Test { void run_cholqr_gaussian(int m, int n, int b, uint32_t key) { DenseDist D(m, n, DenseDistName::Gaussian); std::vector A(m*n); - T iso_scale = RandBLAS::isometry_scale_factor(D); + T iso_scale = D.isometry_scale; RNGState state(key); RandBLAS::fill_dense(D, A.data(), state); blas::scal(m*n, iso_scale, A.data(), 1); @@ -193,7 +193,7 @@ class TestHandrolledQR : public ::testing::Test { void run_qr_blocked_cgs(int m, int n, int b, uint32_t key) { DenseDist D(m, n, DenseDistName::Gaussian); std::vector A(m*n); - T iso_scale = RandBLAS::isometry_scale_factor(D); + T iso_scale = D.isometry_scale; RNGState state(key); RandBLAS::fill_dense(D, A.data(), state); blas::scal(m*n, iso_scale, A.data(), 1); diff --git a/test/test_matmul_cores/test_lskges.cc b/test/test_matmul_cores/test_lskges.cc index 1bbd0d03..25c204f4 100644 --- a/test/test_matmul_cores/test_lskges.cc +++ b/test/test_matmul_cores/test_lskges.cc @@ -57,7 +57,8 @@ class TestLSKGES : public ::testing::Test int64_t nnz_index, int threads ) { - SparseSkOp S0({d, m, vec_nnzs[nnz_index], major_axis}, keys[key_index]); + SparseDist D0(d, m, major_axis, vec_nnzs[nnz_index]); + SparseSkOp S0(D0, keys[key_index]); test_left_apply_to_random(1.0, S0, n, 0.0, layout, threads); } @@ -73,7 +74,8 @@ class TestLSKGES : public ::testing::Test blas::Layout layout ) { int64_t vec_nnz = d0 / 3; // this is actually quite dense. - SparseSkOp S0({d0, m0, vec_nnz, MajorAxis::Short}, seed); + SparseDist D0(d0, m0, MajorAxis::Short, vec_nnz); + SparseSkOp S0(D0, seed); test_left_apply_submatrix_to_eye(1.0, S0, d1, m1, S_ro, S_co, layout, 0.0); } @@ -87,7 +89,7 @@ class TestLSKGES : public ::testing::Test blas::Layout layout ) { int64_t vec_nnz = d / 2; - SparseDist DS = {d, m, vec_nnz, MajorAxis::Short}; + SparseDist DS(d, m, MajorAxis::Short, vec_nnz); SparseSkOp S(DS, key); test_left_apply_submatrix_to_eye(alpha, S, d, m, 0, 0, layout, beta); } @@ -103,12 +105,7 @@ class TestLSKGES : public ::testing::Test randblas_require(m > d); bool is_saso = (major_axis == MajorAxis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; - SparseDist Dt = { - .n_rows = m, - .n_cols = d, - .vec_nnz = vec_nnz, - .major_axis = major_axis - }; + SparseDist Dt(m, d, major_axis, vec_nnz); SparseSkOp S0(Dt, key); test_left_apply_transpose_to_eye(S0, layout); } @@ -129,12 +126,7 @@ class TestLSKGES : public ::testing::Test // Define the distribution for S0. bool is_saso = (major_axis == MajorAxis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; - SparseDist D = { - .n_rows = d, - .n_cols = m, - .vec_nnz = vec_nnz, - .major_axis = major_axis - }; + SparseDist D(d, m, major_axis, vec_nnz); SparseSkOp S0(D, seed_S0); test_left_apply_to_submatrix(S0, n, m0, n0, A_ro, A_co, layout); } @@ -151,12 +143,7 @@ class TestLSKGES : public ::testing::Test // Define the distribution for S0. bool is_saso = (major_axis == MajorAxis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; - SparseDist D = { - .n_rows = d, - .n_cols = m, - .vec_nnz = vec_nnz, - .major_axis = major_axis - }; + SparseDist D(d, m, major_axis, vec_nnz); SparseSkOp S0(D, seed_S0); test_left_apply_to_transposed(S0, n, layout); } diff --git a/test/test_matmul_cores/test_rskges.cc b/test/test_matmul_cores/test_rskges.cc index 35b62f61..7264193a 100644 --- a/test/test_matmul_cores/test_rskges.cc +++ b/test/test_matmul_cores/test_rskges.cc @@ -55,9 +55,7 @@ class TestRSKGES : public ::testing::Test int64_t vec_nnz, Layout layout ) { - SparseDist D = { - .n_rows = m, .n_cols = d, .vec_nnz = vec_nnz, .major_axis = major_axis - }; + SparseDist D(m, d, major_axis, vec_nnz); SparseSkOp S0(D, seed); RandBLAS::fill_sparse(S0); test_right_apply_submatrix_to_eye(1.0, S0, m, d, 0, 0, layout, 0.0, 0); @@ -74,14 +72,10 @@ class TestRSKGES : public ::testing::Test int64_t nnz_index, int threads ) { - SparseDist D = { - .n_rows=n, .n_cols=d, .vec_nnz=vec_nnzs[nnz_index], .major_axis=major_axis - }; + SparseDist D(n, d, major_axis, vec_nnzs[nnz_index]); SparseSkOp S0(D, keys[key_index]); RandBLAS::fill_sparse(S0); test_right_apply_to_random(1.0, S0, m, layout, 0.0, threads); - - } template @@ -98,9 +92,8 @@ class TestRSKGES : public ::testing::Test randblas_require(d0 >= d1); randblas_require(n0 >= n1); int64_t vec_nnz = d0 / 3; // this is actually quite dense. - SparseSkOp S0( - {n0, d0, vec_nnz, RandBLAS::MajorAxis::Short}, seed - ); + SparseDist D0(n0, d0, RandBLAS::MajorAxis::Short, vec_nnz); + SparseSkOp S0(D0, seed); RandBLAS::fill_sparse(S0); test_right_apply_submatrix_to_eye(1.0, S0, n1, d1, S_ro, S_co, layout, 0.0, 0); }