Line data Source code
1 : /** \file
2 : * \brief Reusable row and column deletion updates for thin singular value decompositions.
3 : *
4 : * \ingroup gen_math_files
5 : */
6 :
7 : //***********************************************************************//
8 : // Copyright 2026 Jared R. Males (jaredmales@gmail.com)
9 : //
10 : // This file is part of mxlib.
11 : //
12 : // mxlib is free software: you can redistribute it and/or modify
13 : // it under the terms of the GNU General Public License as published by
14 : // the Free Software Foundation, either version 3 of the License, or
15 : // (at your option) any later version.
16 : //
17 : // mxlib is distributed in the hope that it will be useful,
18 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 : // GNU General Public License for more details.
21 : //
22 : // You should have received a copy of the GNU General Public License
23 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
24 : //***********************************************************************//
25 :
26 : #ifndef math_svdDowndate_hpp
27 : #define math_svdDowndate_hpp
28 :
29 : #include <Eigen/Dense>
30 :
31 : #include <cstddef>
32 : #include <cstdint>
33 : #include <memory>
34 : #include <span>
35 : #include <type_traits>
36 :
37 : #include "templateLapack.hpp"
38 :
39 : namespace mx
40 : {
41 : namespace math
42 : {
43 :
44 : #if defined( __GNUC__ ) || defined( __clang__ )
45 : #define MXLIB_SVD_DELETION_HEADER_ADAPTER inline __attribute__( ( visibility( "hidden" ) ) )
46 : #else
47 : #define MXLIB_SVD_DELETION_HEADER_ADAPTER inline
48 : #endif
49 :
50 : /** \addtogroup svd_downdate
51 : * @{
52 : */
53 :
54 : /// Numerical backend used to delete rows or columns from thin-SVD factors.
55 : enum class svdDeletionBackend
56 : {
57 : leadingCovariance, ///< Symmetric leading-spectrum core; fastest when small singular values are not required.
58 : stableCore, ///< Complement-preserving small SVD; avoids squaring singular-value conditioning.
59 : rankOneSecular ///< Structured covariance eigensolve for deleting exactly one singular-factor row.
60 : };
61 :
62 : /// Completion status for an SVD deletion operation.
63 : enum class svdDeletionStatus
64 : {
65 : notComputed, ///< No operation has published a result.
66 : success, ///< The operation completed without numerical clamping.
67 : successWithClamping, ///< The operation completed after clamping roundoff-scale negative eigenvalues.
68 : invalidInput, ///< Dimensions, values, indices, or requested output rank are invalid.
69 : allocationFailure, ///< Result or workspace allocation failed.
70 : workspaceQueryFailure, ///< LAPACK returned an invalid or failed workspace query.
71 : solverFailure, ///< LAPACK failed during the numerical solve.
72 : nonFiniteOutput, ///< LAPACK returned a non-finite singular system.
73 : invalidSolverOutput, ///< LAPACK returned a finite spectrum with invalid ordering or sign.
74 : rescalingOverflow, ///< A finite normalized result cannot be represented after restoring input scale.
75 : nonPositiveSemidefinite, ///< A theoretically PSD core has a materially negative eigenvalue.
76 : factorNotOrthonormal, ///< A requested singular-factor validation failed.
77 : unsupportedDeletionCount ///< The selected backend cannot process the requested number of deleted rows.
78 : };
79 :
80 : /// Return a stable text representation of an SVD deletion backend.
81 : const char *svdDeletionBackendName( svdDeletionBackend backend /**< [in] backend to describe */ );
82 :
83 : /// Return a stable text representation of an SVD deletion status.
84 : const char *svdDeletionStatusName( svdDeletionStatus status /**< [in] status to describe */ );
85 :
86 : /// Return true when a status represents usable numerical output.
87 : bool svdDeletionSucceeded( svdDeletionStatus status /**< [in] status to classify */ ) noexcept;
88 :
89 : /// Column-major dynamic matrix used by the SVD deletion API.
90 : template <typename realT>
91 : using svdDeletionMatrix = Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
92 :
93 : /// Dynamic column vector used by the SVD deletion API.
94 : template <typename realT>
95 : using svdDeletionVector = Eigen::Array<realT, Eigen::Dynamic, 1>;
96 :
97 : /// Non-owning read-only reference to a compatible column-major SVD deletion matrix.
98 : template <typename realT>
99 : using svdDeletionConstMatrixRef = Eigen::Ref<const svdDeletionMatrix<realT>>;
100 :
101 : /// Non-owning read-only reference to a compatible SVD deletion vector.
102 : template <typename realT>
103 : using svdDeletionConstVectorRef = Eigen::Ref<const svdDeletionVector<realT>>;
104 :
105 : static_assert( std::is_integral_v<Eigen::Index>, "The SVD-deletion ABI requires Eigen::Index to be an integral type." );
106 : static_assert( std::is_signed_v<Eigen::Index>, "The SVD-deletion ABI requires Eigen::Index to be signed." );
107 : static_assert( sizeof( Eigen::Index ) <= sizeof( std::int64_t ),
108 : "The SVD-deletion ABI cannot represent an Eigen::Index wider than int64_t." );
109 : static_assert( sizeof( Eigen::Index ) == sizeof( std::ptrdiff_t ),
110 : "The SVD-deletion ABI requires Eigen::Index and ptrdiff_t to have the same width." );
111 :
112 : /// ABI-stable borrowed contiguous-vector storage descriptor.
113 : template <typename realT>
114 : struct svdDeletionConstVectorViewV2
115 : {
116 : const realT *data{ nullptr }; ///< First scalar, or null for an empty view.
117 :
118 : std::int64_t size{ 0 }; ///< Number of contiguous scalars.
119 : };
120 :
121 : /// ABI-stable borrowed column-major matrix storage descriptor.
122 : template <typename realT>
123 : struct svdDeletionConstMatrixViewV2
124 : {
125 : const realT *data{ nullptr }; ///< First scalar, or null for an empty view.
126 :
127 : std::int64_t rows{ 0 }; ///< Matrix row count.
128 :
129 : std::int64_t columns{ 0 }; ///< Matrix column count.
130 :
131 : std::int64_t outerStride{ 0 }; ///< Scalar stride between successive columns.
132 : };
133 :
134 : /// ABI-stable borrowed signed-index storage descriptor.
135 : struct svdDeletionConstIndexViewV2
136 : {
137 : const void *data{ nullptr }; ///< First signed index, or null for an empty view.
138 :
139 : std::int64_t size{ 0 }; ///< Number of indices.
140 :
141 : std::int64_t elementBytes{ 0 }; ///< Width of each signed integer element.
142 : };
143 :
144 : static_assert( std::is_standard_layout_v<svdDeletionConstVectorViewV2<float>> );
145 : static_assert( std::is_trivially_copyable_v<svdDeletionConstVectorViewV2<float>> );
146 : static_assert( std::is_standard_layout_v<svdDeletionConstVectorViewV2<double>> );
147 : static_assert( std::is_trivially_copyable_v<svdDeletionConstVectorViewV2<double>> );
148 : static_assert( std::is_standard_layout_v<svdDeletionConstMatrixViewV2<float>> );
149 : static_assert( std::is_trivially_copyable_v<svdDeletionConstMatrixViewV2<float>> );
150 : static_assert( std::is_standard_layout_v<svdDeletionConstMatrixViewV2<double>> );
151 : static_assert( std::is_trivially_copyable_v<svdDeletionConstMatrixViewV2<double>> );
152 : static_assert( std::is_standard_layout_v<svdDeletionConstIndexViewV2> );
153 : static_assert( std::is_trivially_copyable_v<svdDeletionConstIndexViewV2> );
154 :
155 : /// Type-level ABI tag for the second-generation opaque SVD deletion handles.
156 : struct svdDeletionAbiV2Tag
157 : {
158 : };
159 :
160 : template <typename realT, typename abiT = svdDeletionAbiV2Tag>
161 : class svdDeletionResult;
162 :
163 : template <typename realT, typename abiT = svdDeletionAbiV2Tag>
164 : class svdDeletionWorkspace;
165 :
166 : namespace detail
167 : {
168 :
169 : template <typename realT>
170 : struct svdDeletionImplementation;
171 :
172 : /// \cond svdDeletion_abi_detail
173 :
174 : // ABI-v2 entry point for factor validation without Eigen types at the shared-library boundary.
175 : svdDeletionStatus validateSvdDeletionFactorAbiV2( svdDeletionConstMatrixViewV2<float> factor, float tolerance );
176 :
177 : // ABI-v2 entry point for double-precision factor validation.
178 : svdDeletionStatus validateSvdDeletionFactorAbiV2( svdDeletionConstMatrixViewV2<double> factor, double tolerance );
179 :
180 : // ABI-v2 entry point for the leading-covariance deletion core.
181 : template <typename realT>
182 : svdDeletionStatus svdDeletionLeadingCoreAbiV2( svdDeletionResult<realT> &result,
183 : svdDeletionConstVectorViewV2<realT> singularValues,
184 : svdDeletionConstMatrixViewV2<realT> deletedRows,
185 : std::int64_t outputRank,
186 : svdDeletionWorkspace<realT> &workspace );
187 :
188 : // ABI-v2 entry point for the stable deletion core.
189 : template <typename realT>
190 : svdDeletionStatus svdDeletionStableCoreAbiV2( svdDeletionResult<realT> &result,
191 : svdDeletionConstVectorViewV2<realT> singularValues,
192 : svdDeletionConstMatrixViewV2<realT> deletedRows,
193 : std::int64_t outputRank,
194 : svdDeletionWorkspace<realT> &workspace );
195 :
196 : // ABI-v2 entry point for backend-selected factor deletion.
197 : template <typename realT>
198 : svdDeletionStatus svdDeletionCoreAbiV2( svdDeletionResult<realT> &result,
199 : svdDeletionConstVectorViewV2<realT> singularValues,
200 : svdDeletionConstMatrixViewV2<realT> deletedRows,
201 : std::int64_t outputRank,
202 : svdDeletionWorkspace<realT> &workspace,
203 : svdDeletionBackend backend );
204 :
205 : // ABI-v2 entry point for physical row deletion.
206 : template <typename realT>
207 : svdDeletionStatus svdRemoveRowsAbiV2( svdDeletionResult<realT> &result,
208 : svdDeletionConstVectorViewV2<realT> singularValues,
209 : svdDeletionConstMatrixViewV2<realT> leftFactor,
210 : svdDeletionConstIndexViewV2 deletedIndices,
211 : std::int64_t outputRank,
212 : svdDeletionWorkspace<realT> &workspace,
213 : svdDeletionBackend backend );
214 :
215 : // ABI-v2 entry point for physical column deletion.
216 : template <typename realT>
217 : svdDeletionStatus svdRemoveColumnsAbiV2( svdDeletionResult<realT> &result,
218 : svdDeletionConstVectorViewV2<realT> singularValues,
219 : svdDeletionConstMatrixViewV2<realT> rightFactor,
220 : svdDeletionConstIndexViewV2 deletedIndices,
221 : std::int64_t outputRank,
222 : svdDeletionWorkspace<realT> &workspace,
223 : svdDeletionBackend backend );
224 :
225 : /// \endcond
226 :
227 : /// \cond svdDeletion_test_detail
228 :
229 : /// Internal stages exposed only for deterministic failure-path tests.
230 : enum class svdDeletionTestOperation
231 : {
232 : validateFactor,
233 : prepareResult,
234 : prepareWorkspace
235 : };
236 :
237 : /// Function signature for deterministic allocation-failure injection.
238 : using svdDeletionOperationHookT = void ( * )( svdDeletionTestOperation );
239 :
240 : /// Function signature matching the LAPACK SYEVR wrapper.
241 : template <typename realT>
242 : using svdDeletionSyevrHookT = MXLAPACK_INT ( * )( char,
243 : char,
244 : char,
245 : MXLAPACK_INT,
246 : realT *,
247 : MXLAPACK_INT,
248 : realT,
249 : realT,
250 : MXLAPACK_INT,
251 : MXLAPACK_INT,
252 : realT,
253 : MXLAPACK_INT *,
254 : realT *,
255 : realT *,
256 : MXLAPACK_INT,
257 : MXLAPACK_INT *,
258 : realT *,
259 : MXLAPACK_INT,
260 : MXLAPACK_INT *,
261 : MXLAPACK_INT );
262 :
263 : /// Function signature matching the LAPACK GESVD wrapper.
264 : template <typename realT>
265 : using svdDeletionGesvdHookT = MXLAPACK_INT ( * )( char,
266 : char,
267 : MXLAPACK_INT,
268 : MXLAPACK_INT,
269 : realT *,
270 : MXLAPACK_INT,
271 : realT *,
272 : realT *,
273 : MXLAPACK_INT,
274 : realT *,
275 : MXLAPACK_INT,
276 : realT *,
277 : MXLAPACK_INT );
278 :
279 : /// Function signature for the structured rank-one secular eigensolver hook.
280 : template <typename realT>
281 : using svdDeletionLaed9HookT =
282 : MXLAPACK_INT ( * )( realT *, realT *, realT *, MXLAPACK_INT, MXLAPACK_INT, realT, realT *, realT * );
283 :
284 : /// Test hooks for one floating-point specialization.
285 : template <typename realT>
286 : struct svdDeletionTestHooks
287 : {
288 : /// Optional operation hook; null performs the production operation.
289 : svdDeletionOperationHookT operation{ nullptr };
290 :
291 : /// Optional SYEVR hook; null calls `math::syevr<realT>`.
292 : svdDeletionSyevrHookT<realT> syevr{ nullptr };
293 :
294 : /// Optional GESVD hook; null calls `math::gesvd<realT>`.
295 : svdDeletionGesvdHookT<realT> gesvd{ nullptr };
296 :
297 : /// Optional LAED9 hook; null calls the LAPACK structured secular eigensolver.
298 : svdDeletionLaed9HookT<realT> laed9{ nullptr };
299 : };
300 :
301 : /// Access the process-wide SVD deletion test hooks for one scalar type.
302 : template <typename realT>
303 : svdDeletionTestHooks<realT> &svdDeletionHooks();
304 :
305 : /// \endcond
306 :
307 : } // namespace detail
308 :
309 : /// Result of deleting rows or columns from a represented thin SVD.
310 : /** If \f$A=U\Sigma V^T\f$, the returned rotation applies to the factor on the side that was not deleted. For row
311 : * deletion the updated right factor is `V * rotation()`; for column deletion the updated left factor is
312 : * `U * rotation()`. Singular values are descending.
313 : *
314 : * Output arrays remain allocated across calls. Their numerical contents are valid only when status() reports a
315 : * successful operation.
316 : *
317 : * \tparam realT floating-point type; supported explicit instantiations are float and double.
318 : * \tparam abiT type-level ABI tag; callers use the default.
319 : */
320 : template <typename realT, typename abiT>
321 : class svdDeletionResult
322 : {
323 : static_assert( std::is_same_v<abiT, svdDeletionAbiV2Tag>,
324 : "svdDeletionResult does not support a caller-selected ABI tag." );
325 :
326 : public:
327 : /// Construct an empty result.
328 : svdDeletionResult();
329 :
330 : /// Release result storage using mxlib's Eigen allocation configuration.
331 : ~svdDeletionResult();
332 :
333 : /// Results cannot be copied across the mxlib ABI boundary.
334 : svdDeletionResult( const svdDeletionResult &other /**< [in] result that copying is forbidden from */ ) = delete;
335 :
336 : /// Results cannot be copy-assigned across the mxlib ABI boundary.
337 : svdDeletionResult &
338 : operator=( const svdDeletionResult &other /**< [in] result that copy assignment is forbidden from */ ) = delete;
339 :
340 : /// Move owned result storage from another result.
341 : svdDeletionResult( svdDeletionResult &&other /**< [in,out] result to move from */ ) noexcept;
342 :
343 : /// Replace this result by moving owned storage from another result.
344 : svdDeletionResult &operator=( svdDeletionResult &&other /**< [in,out] result to move from */ ) noexcept;
345 :
346 : /// Prepare or reuse output storage for a base and requested active output rank.
347 : /** Storage grows when required but does not shrink while `baseRank` is unchanged. Calling this before a hot loop
348 : * with the largest planned output rank reserves the rotation capacity for later smaller requests.
349 : */
350 : MXLIB_SVD_DELETION_HEADER_ADAPTER
351 120 : svdDeletionStatus prepare( Eigen::Index baseRank, /**< [in] number of supplied singular triplets */
352 : Eigen::Index outputRank /**< [in] number of updated triplets to publish */ )
353 : {
354 120 : return prepareAbiV2( static_cast<std::int64_t>( baseRank ), static_cast<std::int64_t>( outputRank ) );
355 : }
356 :
357 : /// Return the most recent operation status.
358 : svdDeletionStatus status() const noexcept;
359 :
360 : /// Return the backend that produced the current result.
361 : svdDeletionBackend backend() const noexcept;
362 :
363 : /// Return an unaligned borrowed view of all `baseRank()` descending updated singular values.
364 : /** The view remains valid until this result is prepared, assigned, moved, or destroyed. */
365 81 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionConstVectorRef<realT> singularValues() const noexcept
366 : {
367 81 : const svdDeletionConstVectorViewV2<realT> view = singularValuesViewAbiV2();
368 162 : return Eigen::Map<const svdDeletionVector<realT>, Eigen::Unaligned>( view.data,
369 81 : static_cast<Eigen::Index>( view.size ) );
370 : }
371 :
372 : /// Return an unaligned borrowed view of all corresponding descending squared singular values.
373 : /** The view remains valid until this result is prepared, assigned, moved, or destroyed. */
374 102 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionConstVectorRef<realT> squaredSingularValues() const noexcept
375 : {
376 102 : const svdDeletionConstVectorViewV2<realT> view = squaredSingularValuesViewAbiV2();
377 204 : return Eigen::Map<const svdDeletionVector<realT>, Eigen::Unaligned>( view.data,
378 102 : static_cast<Eigen::Index>( view.size ) );
379 : }
380 :
381 : /// Return the preserved-side rotation, with updated directions in columns.
382 : /** The unaligned borrowed view remains valid until this result is prepared, assigned, moved, or destroyed. */
383 50 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionConstMatrixRef<realT> rotation() const noexcept
384 : {
385 50 : const svdDeletionConstMatrixViewV2<realT> view = rotationViewAbiV2();
386 : using mapT = Eigen::Map<const svdDeletionMatrix<realT>, Eigen::Unaligned, Eigen::OuterStride<Eigen::Dynamic>>;
387 100 : return mapT( view.data,
388 50 : static_cast<Eigen::Index>( view.rows ),
389 50 : static_cast<Eigen::Index>( view.columns ),
390 100 : Eigen::OuterStride<Eigen::Dynamic>( static_cast<Eigen::Index>( view.outerStride ) ) );
391 : }
392 :
393 : /// Return the base factorization rank for which storage is prepared.
394 : std::int64_t baseRank() const noexcept;
395 :
396 : /// Return the requested published rank.
397 : std::int64_t outputRank() const noexcept;
398 :
399 : /// Return the allocated rotation-column capacity for the current base rank.
400 : std::int64_t maximumOutputRank() const noexcept;
401 :
402 : /// Return the number of roundoff-scale negative eigenvalues clamped to zero.
403 : std::int64_t clampedEigenvalues() const noexcept;
404 :
405 : /// Return the smallest pre-clamp eigenvalue from the normalized backend PSD validation core.
406 : /** A non-positive-semidefinite failure preserves the offending value. An empty stable-core deletion uses one as
407 : * the vacuous minimum of its zero-dimensional complement core.
408 : */
409 : realT minimumPSDValue() const noexcept;
410 :
411 : /// Return the underlying LAPACK status from the most recent failed query or solve.
412 : MXLAPACK_INT lapackInfo() const noexcept;
413 :
414 : private:
415 : friend struct detail::svdDeletionImplementation<realT>;
416 :
417 : /// Prepare output storage through the ABI-v2 integer-only boundary.
418 : svdDeletionStatus prepareAbiV2( std::int64_t baseRank, /**< [in] number of supplied singular triplets */
419 : std::int64_t outputRank /**< [in] number of triplets to publish */ );
420 :
421 : /// Return the singular-value storage through the ABI-v2 POD boundary.
422 : svdDeletionConstVectorViewV2<realT> singularValuesViewAbiV2() const noexcept;
423 :
424 : /// Return squared-singular-value storage through the ABI-v2 POD boundary.
425 : svdDeletionConstVectorViewV2<realT> squaredSingularValuesViewAbiV2() const noexcept;
426 :
427 : /// Return rotation storage through the ABI-v2 POD boundary.
428 : svdDeletionConstMatrixViewV2<realT> rotationViewAbiV2() const noexcept;
429 :
430 : /// Allocate private result storage after a move, translating allocation failure to false.
431 : bool ensureStorage() noexcept;
432 :
433 : class storage;
434 :
435 : /// Opaque result storage allocated and released with mxlib's Eigen build configuration.
436 : std::unique_ptr<storage> m_storage;
437 : };
438 :
439 : /// Reusable, non-shared storage for SVD deletion operations.
440 : /** Call prepare() before a hot loop to allocate and query LAPACK once. A later operation reuses the storage whenever
441 : * its dimensions and backend fit the prepared capacity. The workspace is noncopyable but movable, allowing callers
442 : * to own one independent workspace per worker without sharing in-flight numerical state.
443 : *
444 : * \tparam realT floating-point type; supported explicit instantiations are float and double.
445 : * \tparam abiT type-level ABI tag; callers use the default.
446 : */
447 : template <typename realT, typename abiT>
448 : class svdDeletionWorkspace
449 : {
450 : static_assert( std::is_same_v<abiT, svdDeletionAbiV2Tag>,
451 : "svdDeletionWorkspace does not support a caller-selected ABI tag." );
452 :
453 : public:
454 : /// Construct an empty workspace.
455 : svdDeletionWorkspace();
456 :
457 : /// Release owned scratch storage.
458 : ~svdDeletionWorkspace();
459 :
460 : /// Workspaces cannot be copied.
461 : svdDeletionWorkspace( const svdDeletionWorkspace &other /**< [in] workspace that copying is forbidden from */ ) =
462 : delete;
463 :
464 : /// Workspaces cannot be copy-assigned.
465 : svdDeletionWorkspace &operator=(
466 : const svdDeletionWorkspace &other /**< [in] workspace that copy assignment is forbidden from */ ) = delete;
467 :
468 : /// Move owned storage and preparation state from another workspace.
469 : svdDeletionWorkspace( svdDeletionWorkspace &&other /**< [in,out] workspace to move from */ ) noexcept;
470 :
471 : /// Replace this workspace by moving owned storage and preparation state.
472 : svdDeletionWorkspace &operator=( svdDeletionWorkspace &&other /**< [in,out] workspace to move from */ ) noexcept;
473 :
474 : /// Prepare reusable storage and LAPACK work arrays.
475 : MXLIB_SVD_DELETION_HEADER_ADAPTER
476 103 : svdDeletionStatus prepare( Eigen::Index baseRank, /**< [in] maximum supplied singular rank */
477 : Eigen::Index maximumDeleted, /**< [in] maximum rows of the deleted-side factor */
478 : svdDeletionBackend backend /**< [in] numerical backend to prepare */ )
479 : {
480 103 : return prepareAbiV2( static_cast<std::int64_t>( baseRank ),
481 : static_cast<std::int64_t>( maximumDeleted ),
482 103 : backend );
483 : }
484 :
485 : /// Release all prepared storage and reset dimensions.
486 : void clear() noexcept;
487 :
488 : /// Report whether this workspace has completed preparation.
489 : bool prepared() const noexcept;
490 :
491 : /// Return the prepared base rank.
492 : std::int64_t baseRank() const noexcept;
493 :
494 : /// Return the prepared maximum deletion count.
495 : std::int64_t maximumDeleted() const noexcept;
496 :
497 : /// Return the prepared numerical backend.
498 : svdDeletionBackend backend() const noexcept;
499 :
500 : /// Return the underlying LAPACK status from the most recent failed workspace query.
501 : MXLAPACK_INT lapackInfo() const noexcept;
502 :
503 : private:
504 : friend struct detail::svdDeletionImplementation<realT>;
505 :
506 : /// Prepare scratch storage through the ABI-v2 integer-only boundary.
507 : svdDeletionStatus prepareAbiV2( std::int64_t baseRank, /**< [in] maximum supplied singular rank */
508 : std::int64_t maximumDeleted, /**< [in] maximum deleted-factor rows */
509 : svdDeletionBackend backend /**< [in] numerical backend */ );
510 :
511 : /// Allocate private workspace storage after a move, translating allocation failure to false.
512 : bool ensureStorage() noexcept;
513 :
514 : class storage;
515 :
516 : /// Opaque scratch storage allocated and released with mxlib's Eigen build configuration.
517 : std::unique_ptr<storage> m_storage;
518 : };
519 :
520 : /// Validate that a supplied thin singular-vector factor has orthonormal columns.
521 : /** This is an optional one-time base-factor check. Hot deletion calls assume the factor contract and do not repeat
522 : * this `O(n q^2)` operation. A zero tolerance selects a dimension-scaled default.
523 : */
524 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus
525 1 : validateSvdDeletionFactor( svdDeletionConstMatrixRef<float> factor, /**< [in] thin singular-vector factor to validate */
526 : float tolerance = 0 /**< [in] maximum absolute Gram-matrix error, or zero for automatic */ )
527 : {
528 1 : return detail::validateSvdDeletionFactorAbiV2( { factor.data(),
529 1 : static_cast<std::int64_t>( factor.rows() ),
530 1 : static_cast<std::int64_t>( factor.cols() ),
531 1 : static_cast<std::int64_t>( factor.outerStride() ) },
532 1 : tolerance );
533 : }
534 :
535 : /// Validate that a supplied double-precision thin singular-vector factor has orthonormal columns.
536 13 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus validateSvdDeletionFactor(
537 : svdDeletionConstMatrixRef<double> factor, /**< [in] thin singular-vector factor to validate */
538 : double tolerance = 0 /**< [in] maximum absolute Gram-matrix error, or zero for automatic */ )
539 : {
540 13 : return detail::validateSvdDeletionFactorAbiV2( { factor.data(),
541 13 : static_cast<std::int64_t>( factor.rows() ),
542 13 : static_cast<std::int64_t>( factor.cols() ),
543 13 : static_cast<std::int64_t>( factor.outerStride() ) },
544 13 : tolerance );
545 : }
546 :
547 : /// Delete supplied singular-factor rows with the full-spectrum symmetric covariance core.
548 : /** Given deleted-side rows `F`, this solves
549 : *
550 : * \f[
551 : * H = \Sigma (I-F^T F) \Sigma = W \Lambda W^T.
552 : * \f]
553 : *
554 : * The result rotation is `W` and its singular values are `sqrt(diag(Lambda))`. The complete spectrum is evaluated
555 : * for PSD validation even when only `outputRank` leading directions are published. This backend forms normal
556 : * equations and therefore does not promise high relative accuracy for the smallest singular values.
557 : * The `rankOneSecular` backend evaluates this same core in quadratic time when exactly one row is deleted.
558 : */
559 : template <typename realT>
560 12 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus svdDeletionLeadingCore(
561 : svdDeletionResult<realT> &result, /**< [out] updated spectrum, rotation, and diagnostics */
562 : std::type_identity_t<svdDeletionConstVectorRef<realT>> singularValues, /**< [in] descending base singular values */
563 : std::type_identity_t<svdDeletionConstMatrixRef<realT>> deletedRows, /**< [in] deleted-side factor rows */
564 : Eigen::Index outputRank, /**< [in] number of leading updated directions to publish */
565 : svdDeletionWorkspace<realT> &workspace /**< [in,out] reusable, worker-private scratch storage */ )
566 : {
567 12 : return detail::svdDeletionLeadingCoreAbiV2<realT>(
568 : result,
569 12 : { singularValues.data(), static_cast<std::int64_t>( singularValues.size() ) },
570 12 : { deletedRows.data(),
571 12 : static_cast<std::int64_t>( deletedRows.rows() ),
572 12 : static_cast<std::int64_t>( deletedRows.cols() ),
573 12 : static_cast<std::int64_t>( deletedRows.outerStride() ) },
574 : static_cast<std::int64_t>( outputRank ),
575 12 : workspace );
576 : }
577 :
578 : /// Delete supplied singular-factor rows with the complement-preserving small-SVD core.
579 : /** Let `F` contain deleted-side singular-factor rows and choose `B` so that
580 : * `B^T B = I - F F^T`. The at-most `(q+c) x q` core
581 : *
582 : * \f[
583 : * K = \begin{bmatrix}(I-F^TF)\Sigma \\ BF\Sigma\end{bmatrix}
584 : * \f]
585 : *
586 : * obeys `K^T K = Sigma (I-F^T F) Sigma`. Its right singular vectors are the preserved-side rotation. This is the
587 : * default generic backend because it avoids explicitly squaring the represented singular values. See
588 : * \cite brand_2006 and \cite long_males_2021.
589 : *
590 : * This backend and `leadingCovariance` use dense LAPACK solvers and therefore have cubic asymptotic cost in the base
591 : * rank. The `rankOneSecular` backend instead uses a structured quadratic-time solve for one-row deletion.
592 : */
593 : template <typename realT>
594 14 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus svdDeletionStableCore(
595 : svdDeletionResult<realT> &result, /**< [out] updated spectrum, rotation, and diagnostics */
596 : std::type_identity_t<svdDeletionConstVectorRef<realT>> singularValues, /**< [in] descending base singular values */
597 : std::type_identity_t<svdDeletionConstMatrixRef<realT>> deletedRows, /**< [in] deleted-side factor rows */
598 : Eigen::Index outputRank, /**< [in] number of leading updated directions to publish */
599 : svdDeletionWorkspace<realT> &workspace /**< [in,out] reusable, worker-private scratch storage */ )
600 : {
601 14 : return detail::svdDeletionStableCoreAbiV2<realT>(
602 : result,
603 14 : { singularValues.data(), static_cast<std::int64_t>( singularValues.size() ) },
604 14 : { deletedRows.data(),
605 14 : static_cast<std::int64_t>( deletedRows.rows() ),
606 14 : static_cast<std::int64_t>( deletedRows.cols() ),
607 14 : static_cast<std::int64_t>( deletedRows.outerStride() ) },
608 : static_cast<std::int64_t>( outputRank ),
609 14 : workspace );
610 : }
611 :
612 : /// Delete supplied singular-factor rows with an explicitly selected backend.
613 : /** `rankOneSecular` accepts either no deleted rows, which returns the identity update, or exactly one deleted row.
614 : * It solves the same covariance core as `leadingCovariance` by transforming the diagonal-minus-rank-one problem to a
615 : * positive rank-one secular equation. LAPACK-style deflation handles negligible update components and clustered
616 : * poles at a roundoff-scaled tolerance; post-solve validation uses dimension-scaled bounds. See
617 : * \cite bunch_nielsen_1978 and \cite gu_eisenstat_1995.
618 : */
619 : template <typename realT>
620 32 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus svdDeletionCore(
621 : svdDeletionResult<realT> &result, /**< [out] updated spectrum, rotation, and diagnostics */
622 : std::type_identity_t<svdDeletionConstVectorRef<realT>> singularValues, /**< [in] descending base singular values */
623 : std::type_identity_t<svdDeletionConstMatrixRef<realT>> deletedRows, /**< [in] deleted-side factor rows */
624 : Eigen::Index outputRank, /**< [in] number of leading updated directions to publish */
625 : svdDeletionWorkspace<realT> &workspace, /**< [in,out] reusable, worker-private scratch storage */
626 : svdDeletionBackend backend = svdDeletionBackend::stableCore /**< [in] numerical backend */ )
627 : {
628 32 : return detail::svdDeletionCoreAbiV2<realT>(
629 : result,
630 32 : { singularValues.data(), static_cast<std::int64_t>( singularValues.size() ) },
631 32 : { deletedRows.data(),
632 32 : static_cast<std::int64_t>( deletedRows.rows() ),
633 32 : static_cast<std::int64_t>( deletedRows.cols() ),
634 32 : static_cast<std::int64_t>( deletedRows.outerStride() ) },
635 : static_cast<std::int64_t>( outputRank ),
636 : workspace,
637 32 : backend );
638 : }
639 :
640 : /// Delete physical rows from the matrix represented by a thin SVD.
641 : /** For `A=U Sigma V^T`, this gathers `U[deletedIndices,:]`. The returned rotation applies to `V`. Supplying complete
642 : * factors makes the deletion identical to a direct SVD of the physically retained matrix. Supplying truncated
643 : * factors deletes rows exactly from that represented low-rank matrix.
644 : */
645 : template <typename realT>
646 56 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus svdRemoveRows(
647 : svdDeletionResult<realT> &result, /**< [out] updated spectrum, rotation, and diagnostics */
648 : std::type_identity_t<svdDeletionConstVectorRef<realT>> singularValues, /**< [in] descending base singular values */
649 : std::type_identity_t<svdDeletionConstMatrixRef<realT>> leftFactor, /**< [in] thin left factor `U` */
650 : std::span<const Eigen::Index> deletedIndices, /**< [in] sorted, unique row indices to delete */
651 : Eigen::Index outputRank, /**< [in] number of leading updated directions */
652 : svdDeletionWorkspace<realT> &workspace, /**< [in,out] reusable, worker-private scratch */
653 : svdDeletionBackend backend = svdDeletionBackend::stableCore /**< [in] numerical backend */ )
654 : {
655 56 : return detail::svdRemoveRowsAbiV2<realT>(
656 : result,
657 56 : { singularValues.data(), static_cast<std::int64_t>( singularValues.size() ) },
658 56 : { leftFactor.data(),
659 56 : static_cast<std::int64_t>( leftFactor.rows() ),
660 56 : static_cast<std::int64_t>( leftFactor.cols() ),
661 56 : static_cast<std::int64_t>( leftFactor.outerStride() ) },
662 56 : { deletedIndices.data(),
663 56 : static_cast<std::int64_t>( deletedIndices.size() ),
664 : static_cast<std::int64_t>( sizeof( Eigen::Index ) ) },
665 : static_cast<std::int64_t>( outputRank ),
666 : workspace,
667 56 : backend );
668 : }
669 :
670 : /// Delete physical columns from the matrix represented by a thin SVD.
671 : /** For `A=U Sigma V^T`, this gathers `V[deletedIndices,:]`. The returned rotation applies to `U`. Supplying complete
672 : * factors makes the deletion identical to a direct SVD of the physically retained matrix. Supplying truncated
673 : * factors deletes columns exactly from that represented low-rank matrix.
674 : */
675 : template <typename realT>
676 10 : MXLIB_SVD_DELETION_HEADER_ADAPTER svdDeletionStatus svdRemoveColumns(
677 : svdDeletionResult<realT> &result, /**< [out] updated spectrum, rotation, and diagnostics */
678 : std::type_identity_t<svdDeletionConstVectorRef<realT>> singularValues, /**< [in] descending base singular values */
679 : std::type_identity_t<svdDeletionConstMatrixRef<realT>> rightFactor, /**< [in] thin right factor `V` */
680 : std::span<const Eigen::Index> deletedIndices, /**< [in] sorted, unique column indices to delete */
681 : Eigen::Index outputRank, /**< [in] number of leading updated directions */
682 : svdDeletionWorkspace<realT> &workspace, /**< [in,out] reusable, worker-private scratch */
683 : svdDeletionBackend backend = svdDeletionBackend::stableCore /**< [in] numerical backend */ )
684 : {
685 10 : return detail::svdRemoveColumnsAbiV2<realT>(
686 : result,
687 10 : { singularValues.data(), static_cast<std::int64_t>( singularValues.size() ) },
688 10 : { rightFactor.data(),
689 10 : static_cast<std::int64_t>( rightFactor.rows() ),
690 10 : static_cast<std::int64_t>( rightFactor.cols() ),
691 10 : static_cast<std::int64_t>( rightFactor.outerStride() ) },
692 10 : { deletedIndices.data(),
693 10 : static_cast<std::int64_t>( deletedIndices.size() ),
694 : static_cast<std::int64_t>( sizeof( Eigen::Index ) ) },
695 : static_cast<std::int64_t>( outputRank ),
696 : workspace,
697 10 : backend );
698 : }
699 :
700 : /** @} */
701 :
702 : #undef MXLIB_SVD_DELETION_HEADER_ADAPTER
703 :
704 : } // namespace math
705 : } // namespace mx
706 :
707 : #endif // math_svdDowndate_hpp
|
