Coverage for /usr/lib/python3/dist-packages/scipy/linalg/interpolative.py: 11%
188 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#******************************************************************************
2# Copyright (C) 2013 Kenneth L. Ho
3#
4# Redistribution and use in source and binary forms, with or without
5# modification, are permitted provided that the following conditions are met:
6#
7# Redistributions of source code must retain the above copyright notice, this
8# list of conditions and the following disclaimer. Redistributions in binary
9# form must reproduce the above copyright notice, this list of conditions and
10# the following disclaimer in the documentation and/or other materials
11# provided with the distribution.
12#
13# None of the names of the copyright holders may be used to endorse or
14# promote products derived from this software without specific prior written
15# permission.
16#
17# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27# POSSIBILITY OF SUCH DAMAGE.
28#******************************************************************************
30# Python module for interfacing with `id_dist`.
32r"""
33======================================================================
34Interpolative matrix decomposition (:mod:`scipy.linalg.interpolative`)
35======================================================================
37.. moduleauthor:: Kenneth L. Ho <klho@stanford.edu>
39.. versionadded:: 0.13
41.. currentmodule:: scipy.linalg.interpolative
43An interpolative decomposition (ID) of a matrix :math:`A \in
44\mathbb{C}^{m \times n}` of rank :math:`k \leq \min \{ m, n \}` is a
45factorization
47.. math::
48 A \Pi =
49 \begin{bmatrix}
50 A \Pi_{1} & A \Pi_{2}
51 \end{bmatrix} =
52 A \Pi_{1}
53 \begin{bmatrix}
54 I & T
55 \end{bmatrix},
57where :math:`\Pi = [\Pi_{1}, \Pi_{2}]` is a permutation matrix with
58:math:`\Pi_{1} \in \{ 0, 1 \}^{n \times k}`, i.e., :math:`A \Pi_{2} =
59A \Pi_{1} T`. This can equivalently be written as :math:`A = BP`,
60where :math:`B = A \Pi_{1}` and :math:`P = [I, T] \Pi^{\mathsf{T}}`
61are the *skeleton* and *interpolation matrices*, respectively.
63If :math:`A` does not have exact rank :math:`k`, then there exists an
64approximation in the form of an ID such that :math:`A = BP + E`, where
65:math:`\| E \| \sim \sigma_{k + 1}` is on the order of the :math:`(k +
661)`-th largest singular value of :math:`A`. Note that :math:`\sigma_{k
67+ 1}` is the best possible error for a rank-:math:`k` approximation
68and, in fact, is achieved by the singular value decomposition (SVD)
69:math:`A \approx U S V^{*}`, where :math:`U \in \mathbb{C}^{m \times
70k}` and :math:`V \in \mathbb{C}^{n \times k}` have orthonormal columns
71and :math:`S = \mathop{\mathrm{diag}} (\sigma_{i}) \in \mathbb{C}^{k
72\times k}` is diagonal with nonnegative entries. The principal
73advantages of using an ID over an SVD are that:
75- it is cheaper to construct;
76- it preserves the structure of :math:`A`; and
77- it is more efficient to compute with in light of the identity submatrix of :math:`P`.
79Routines
80========
82Main functionality:
84.. autosummary::
85 :toctree: generated/
87 interp_decomp
88 reconstruct_matrix_from_id
89 reconstruct_interp_matrix
90 reconstruct_skel_matrix
91 id_to_svd
92 svd
93 estimate_spectral_norm
94 estimate_spectral_norm_diff
95 estimate_rank
97Support functions:
99.. autosummary::
100 :toctree: generated/
102 seed
103 rand
106References
107==========
109This module uses the ID software package [1]_ by Martinsson, Rokhlin,
110Shkolnisky, and Tygert, which is a Fortran library for computing IDs
111using various algorithms, including the rank-revealing QR approach of
112[2]_ and the more recent randomized methods described in [3]_, [4]_,
113and [5]_. This module exposes its functionality in a way convenient
114for Python users. Note that this module does not add any functionality
115beyond that of organizing a simpler and more consistent interface.
117We advise the user to consult also the `documentation for the ID package
118<http://tygert.com/id_doc.4.pdf>`_.
120.. [1] P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert. "ID: a
121 software package for low-rank approximation of matrices via interpolative
122 decompositions, version 0.2." http://tygert.com/id_doc.4.pdf.
124.. [2] H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. "On the
125 compression of low rank matrices." *SIAM J. Sci. Comput.* 26 (4): 1389--1404,
126 2005. :doi:`10.1137/030602678`.
128.. [3] E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, M.
129 Tygert. "Randomized algorithms for the low-rank approximation of matrices."
130 *Proc. Natl. Acad. Sci. U.S.A.* 104 (51): 20167--20172, 2007.
131 :doi:`10.1073/pnas.0709640104`.
133.. [4] P.G. Martinsson, V. Rokhlin, M. Tygert. "A randomized
134 algorithm for the decomposition of matrices." *Appl. Comput. Harmon. Anal.* 30
135 (1): 47--68, 2011. :doi:`10.1016/j.acha.2010.02.003`.
137.. [5] F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert. "A fast
138 randomized algorithm for the approximation of matrices." *Appl. Comput.
139 Harmon. Anal.* 25 (3): 335--366, 2008. :doi:`10.1016/j.acha.2007.12.002`.
142Tutorial
143========
145Initializing
146------------
148The first step is to import :mod:`scipy.linalg.interpolative` by issuing the
149command:
151>>> import scipy.linalg.interpolative as sli
153Now let's build a matrix. For this, we consider a Hilbert matrix, which is well
154know to have low rank:
156>>> from scipy.linalg import hilbert
157>>> n = 1000
158>>> A = hilbert(n)
160We can also do this explicitly via:
162>>> import numpy as np
163>>> n = 1000
164>>> A = np.empty((n, n), order='F')
165>>> for j in range(n):
166>>> for i in range(n):
167>>> A[i,j] = 1. / (i + j + 1)
169Note the use of the flag ``order='F'`` in :func:`numpy.empty`. This
170instantiates the matrix in Fortran-contiguous order and is important for
171avoiding data copying when passing to the backend.
173We then define multiplication routines for the matrix by regarding it as a
174:class:`scipy.sparse.linalg.LinearOperator`:
176>>> from scipy.sparse.linalg import aslinearoperator
177>>> L = aslinearoperator(A)
179This automatically sets up methods describing the action of the matrix and its
180adjoint on a vector.
182Computing an ID
183---------------
185We have several choices of algorithm to compute an ID. These fall largely
186according to two dichotomies:
1881. how the matrix is represented, i.e., via its entries or via its action on a
189 vector; and
1902. whether to approximate it to a fixed relative precision or to a fixed rank.
192We step through each choice in turn below.
194In all cases, the ID is represented by three parameters:
1961. a rank ``k``;
1972. an index array ``idx``; and
1983. interpolation coefficients ``proj``.
200The ID is specified by the relation
201``np.dot(A[:,idx[:k]], proj) == A[:,idx[k:]]``.
203From matrix entries
204...................
206We first consider a matrix given in terms of its entries.
208To compute an ID to a fixed precision, type:
210>>> k, idx, proj = sli.interp_decomp(A, eps)
212where ``eps < 1`` is the desired precision.
214To compute an ID to a fixed rank, use:
216>>> idx, proj = sli.interp_decomp(A, k)
218where ``k >= 1`` is the desired rank.
220Both algorithms use random sampling and are usually faster than the
221corresponding older, deterministic algorithms, which can be accessed via the
222commands:
224>>> k, idx, proj = sli.interp_decomp(A, eps, rand=False)
226and:
228>>> idx, proj = sli.interp_decomp(A, k, rand=False)
230respectively.
232From matrix action
233..................
235Now consider a matrix given in terms of its action on a vector as a
236:class:`scipy.sparse.linalg.LinearOperator`.
238To compute an ID to a fixed precision, type:
240>>> k, idx, proj = sli.interp_decomp(L, eps)
242To compute an ID to a fixed rank, use:
244>>> idx, proj = sli.interp_decomp(L, k)
246These algorithms are randomized.
248Reconstructing an ID
249--------------------
251The ID routines above do not output the skeleton and interpolation matrices
252explicitly but instead return the relevant information in a more compact (and
253sometimes more useful) form. To build these matrices, write:
255>>> B = sli.reconstruct_skel_matrix(A, k, idx)
257for the skeleton matrix and:
259>>> P = sli.reconstruct_interp_matrix(idx, proj)
261for the interpolation matrix. The ID approximation can then be computed as:
263>>> C = np.dot(B, P)
265This can also be constructed directly using:
267>>> C = sli.reconstruct_matrix_from_id(B, idx, proj)
269without having to first compute ``P``.
271Alternatively, this can be done explicitly as well using:
273>>> B = A[:,idx[:k]]
274>>> P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
275>>> C = np.dot(B, P)
277Computing an SVD
278----------------
280An ID can be converted to an SVD via the command:
282>>> U, S, V = sli.id_to_svd(B, idx, proj)
284The SVD approximation is then:
286>>> C = np.dot(U, np.dot(np.diag(S), np.dot(V.conj().T)))
288The SVD can also be computed "fresh" by combining both the ID and conversion
289steps into one command. Following the various ID algorithms above, there are
290correspondingly various SVD algorithms that one can employ.
292From matrix entries
293...................
295We consider first SVD algorithms for a matrix given in terms of its entries.
297To compute an SVD to a fixed precision, type:
299>>> U, S, V = sli.svd(A, eps)
301To compute an SVD to a fixed rank, use:
303>>> U, S, V = sli.svd(A, k)
305Both algorithms use random sampling; for the determinstic versions, issue the
306keyword ``rand=False`` as above.
308From matrix action
309..................
311Now consider a matrix given in terms of its action on a vector.
313To compute an SVD to a fixed precision, type:
315>>> U, S, V = sli.svd(L, eps)
317To compute an SVD to a fixed rank, use:
319>>> U, S, V = sli.svd(L, k)
321Utility routines
322----------------
324Several utility routines are also available.
326To estimate the spectral norm of a matrix, use:
328>>> snorm = sli.estimate_spectral_norm(A)
330This algorithm is based on the randomized power method and thus requires only
331matrix-vector products. The number of iterations to take can be set using the
332keyword ``its`` (default: ``its=20``). The matrix is interpreted as a
333:class:`scipy.sparse.linalg.LinearOperator`, but it is also valid to supply it
334as a :class:`numpy.ndarray`, in which case it is trivially converted using
335:func:`scipy.sparse.linalg.aslinearoperator`.
337The same algorithm can also estimate the spectral norm of the difference of two
338matrices ``A1`` and ``A2`` as follows:
340>>> diff = sli.estimate_spectral_norm_diff(A1, A2)
342This is often useful for checking the accuracy of a matrix approximation.
344Some routines in :mod:`scipy.linalg.interpolative` require estimating the rank
345of a matrix as well. This can be done with either:
347>>> k = sli.estimate_rank(A, eps)
349or:
351>>> k = sli.estimate_rank(L, eps)
353depending on the representation. The parameter ``eps`` controls the definition
354of the numerical rank.
356Finally, the random number generation required for all randomized routines can
357be controlled via :func:`scipy.linalg.interpolative.seed`. To reset the seed
358values to their original values, use:
360>>> sli.seed('default')
362To specify the seed values, use:
364>>> sli.seed(s)
366where ``s`` must be an integer or array of 55 floats. If an integer, the array
367of floats is obtained by using ``numpy.random.rand`` with the given integer
368seed.
370To simply generate some random numbers, type:
372>>> sli.rand(n)
374where ``n`` is the number of random numbers to generate.
376Remarks
377-------
379The above functions all automatically detect the appropriate interface and work
380with both real and complex data types, passing input arguments to the proper
381backend routine.
383"""
385import scipy.linalg._interpolative_backend as _backend
386import numpy as np
387import sys
389__all__ = [
390 'estimate_rank',
391 'estimate_spectral_norm',
392 'estimate_spectral_norm_diff',
393 'id_to_svd',
394 'interp_decomp',
395 'rand',
396 'reconstruct_interp_matrix',
397 'reconstruct_matrix_from_id',
398 'reconstruct_skel_matrix',
399 'seed',
400 'svd',
401]
403_DTYPE_ERROR = ValueError("invalid input dtype (input must be float64 or complex128)")
404_TYPE_ERROR = TypeError("invalid input type (must be array or LinearOperator)")
405_32BIT_ERROR = ValueError("interpolative decomposition on 32-bit systems "
406 "with complex128 is buggy")
407_IS_32BIT = (sys.maxsize < 2**32)
410def _is_real(A):
411 try:
412 if A.dtype == np.complex128:
413 return False
414 elif A.dtype == np.float64:
415 return True
416 else:
417 raise _DTYPE_ERROR
418 except AttributeError as e:
419 raise _TYPE_ERROR from e
422def seed(seed=None):
423 """
424 Seed the internal random number generator used in this ID package.
426 The generator is a lagged Fibonacci method with 55-element internal state.
428 Parameters
429 ----------
430 seed : int, sequence, 'default', optional
431 If 'default', the random seed is reset to a default value.
433 If `seed` is a sequence containing 55 floating-point numbers
434 in range [0,1], these are used to set the internal state of
435 the generator.
437 If the value is an integer, the internal state is obtained
438 from `numpy.random.RandomState` (MT19937) with the integer
439 used as the initial seed.
441 If `seed` is omitted (None), ``numpy.random.rand`` is used to
442 initialize the generator.
444 """
445 # For details, see :func:`_backend.id_srand`, :func:`_backend.id_srandi`,
446 # and :func:`_backend.id_srando`.
448 if isinstance(seed, str) and seed == 'default':
449 _backend.id_srando()
450 elif hasattr(seed, '__len__'):
451 state = np.asfortranarray(seed, dtype=float)
452 if state.shape != (55,):
453 raise ValueError("invalid input size")
454 elif state.min() < 0 or state.max() > 1:
455 raise ValueError("values not in range [0,1]")
456 _backend.id_srandi(state)
457 elif seed is None:
458 _backend.id_srandi(np.random.rand(55))
459 else:
460 rnd = np.random.RandomState(seed)
461 _backend.id_srandi(rnd.rand(55))
464def rand(*shape):
465 """
466 Generate standard uniform pseudorandom numbers via a very efficient lagged
467 Fibonacci method.
469 This routine is used for all random number generation in this package and
470 can affect ID and SVD results.
472 Parameters
473 ----------
474 *shape
475 Shape of output array
477 """
478 # For details, see :func:`_backend.id_srand`, and :func:`_backend.id_srando`.
479 return _backend.id_srand(np.prod(shape)).reshape(shape)
482def interp_decomp(A, eps_or_k, rand=True):
483 """
484 Compute ID of a matrix.
486 An ID of a matrix `A` is a factorization defined by a rank `k`, a column
487 index array `idx`, and interpolation coefficients `proj` such that::
489 numpy.dot(A[:,idx[:k]], proj) = A[:,idx[k:]]
491 The original matrix can then be reconstructed as::
493 numpy.hstack([A[:,idx[:k]],
494 numpy.dot(A[:,idx[:k]], proj)]
495 )[:,numpy.argsort(idx)]
497 or via the routine :func:`reconstruct_matrix_from_id`. This can
498 equivalently be written as::
500 numpy.dot(A[:,idx[:k]],
501 numpy.hstack([numpy.eye(k), proj])
502 )[:,np.argsort(idx)]
504 in terms of the skeleton and interpolation matrices::
506 B = A[:,idx[:k]]
508 and::
510 P = numpy.hstack([numpy.eye(k), proj])[:,np.argsort(idx)]
512 respectively. See also :func:`reconstruct_interp_matrix` and
513 :func:`reconstruct_skel_matrix`.
515 The ID can be computed to any relative precision or rank (depending on the
516 value of `eps_or_k`). If a precision is specified (`eps_or_k < 1`), then
517 this function has the output signature::
519 k, idx, proj = interp_decomp(A, eps_or_k)
521 Otherwise, if a rank is specified (`eps_or_k >= 1`), then the output
522 signature is::
524 idx, proj = interp_decomp(A, eps_or_k)
526 .. This function automatically detects the form of the input parameters
527 and passes them to the appropriate backend. For details, see
528 :func:`_backend.iddp_id`, :func:`_backend.iddp_aid`,
529 :func:`_backend.iddp_rid`, :func:`_backend.iddr_id`,
530 :func:`_backend.iddr_aid`, :func:`_backend.iddr_rid`,
531 :func:`_backend.idzp_id`, :func:`_backend.idzp_aid`,
532 :func:`_backend.idzp_rid`, :func:`_backend.idzr_id`,
533 :func:`_backend.idzr_aid`, and :func:`_backend.idzr_rid`.
535 Parameters
536 ----------
537 A : :class:`numpy.ndarray` or :class:`scipy.sparse.linalg.LinearOperator` with `rmatvec`
538 Matrix to be factored
539 eps_or_k : float or int
540 Relative error (if `eps_or_k < 1`) or rank (if `eps_or_k >= 1`) of
541 approximation.
542 rand : bool, optional
543 Whether to use random sampling if `A` is of type :class:`numpy.ndarray`
544 (randomized algorithms are always used if `A` is of type
545 :class:`scipy.sparse.linalg.LinearOperator`).
547 Returns
548 -------
549 k : int
550 Rank required to achieve specified relative precision if
551 `eps_or_k < 1`.
552 idx : :class:`numpy.ndarray`
553 Column index array.
554 proj : :class:`numpy.ndarray`
555 Interpolation coefficients.
556 """
557 from scipy.sparse.linalg import LinearOperator
559 real = _is_real(A)
561 if isinstance(A, np.ndarray):
562 if eps_or_k < 1:
563 eps = eps_or_k
564 if rand:
565 if real:
566 k, idx, proj = _backend.iddp_aid(eps, A)
567 else:
568 if _IS_32BIT:
569 raise _32BIT_ERROR
570 k, idx, proj = _backend.idzp_aid(eps, A)
571 else:
572 if real:
573 k, idx, proj = _backend.iddp_id(eps, A)
574 else:
575 k, idx, proj = _backend.idzp_id(eps, A)
576 return k, idx - 1, proj
577 else:
578 k = int(eps_or_k)
579 if rand:
580 if real:
581 idx, proj = _backend.iddr_aid(A, k)
582 else:
583 if _IS_32BIT:
584 raise _32BIT_ERROR
585 idx, proj = _backend.idzr_aid(A, k)
586 else:
587 if real:
588 idx, proj = _backend.iddr_id(A, k)
589 else:
590 idx, proj = _backend.idzr_id(A, k)
591 return idx - 1, proj
592 elif isinstance(A, LinearOperator):
593 m, n = A.shape
594 matveca = A.rmatvec
595 if eps_or_k < 1:
596 eps = eps_or_k
597 if real:
598 k, idx, proj = _backend.iddp_rid(eps, m, n, matveca)
599 else:
600 if _IS_32BIT:
601 raise _32BIT_ERROR
602 k, idx, proj = _backend.idzp_rid(eps, m, n, matveca)
603 return k, idx - 1, proj
604 else:
605 k = int(eps_or_k)
606 if real:
607 idx, proj = _backend.iddr_rid(m, n, matveca, k)
608 else:
609 if _IS_32BIT:
610 raise _32BIT_ERROR
611 idx, proj = _backend.idzr_rid(m, n, matveca, k)
612 return idx - 1, proj
613 else:
614 raise _TYPE_ERROR
617def reconstruct_matrix_from_id(B, idx, proj):
618 """
619 Reconstruct matrix from its ID.
621 A matrix `A` with skeleton matrix `B` and ID indices and coefficients `idx`
622 and `proj`, respectively, can be reconstructed as::
624 numpy.hstack([B, numpy.dot(B, proj)])[:,numpy.argsort(idx)]
626 See also :func:`reconstruct_interp_matrix` and
627 :func:`reconstruct_skel_matrix`.
629 .. This function automatically detects the matrix data type and calls the
630 appropriate backend. For details, see :func:`_backend.idd_reconid` and
631 :func:`_backend.idz_reconid`.
633 Parameters
634 ----------
635 B : :class:`numpy.ndarray`
636 Skeleton matrix.
637 idx : :class:`numpy.ndarray`
638 Column index array.
639 proj : :class:`numpy.ndarray`
640 Interpolation coefficients.
642 Returns
643 -------
644 :class:`numpy.ndarray`
645 Reconstructed matrix.
646 """
647 if _is_real(B):
648 return _backend.idd_reconid(B, idx + 1, proj)
649 else:
650 return _backend.idz_reconid(B, idx + 1, proj)
653def reconstruct_interp_matrix(idx, proj):
654 """
655 Reconstruct interpolation matrix from ID.
657 The interpolation matrix can be reconstructed from the ID indices and
658 coefficients `idx` and `proj`, respectively, as::
660 P = numpy.hstack([numpy.eye(proj.shape[0]), proj])[:,numpy.argsort(idx)]
662 The original matrix can then be reconstructed from its skeleton matrix `B`
663 via::
665 numpy.dot(B, P)
667 See also :func:`reconstruct_matrix_from_id` and
668 :func:`reconstruct_skel_matrix`.
670 .. This function automatically detects the matrix data type and calls the
671 appropriate backend. For details, see :func:`_backend.idd_reconint` and
672 :func:`_backend.idz_reconint`.
674 Parameters
675 ----------
676 idx : :class:`numpy.ndarray`
677 Column index array.
678 proj : :class:`numpy.ndarray`
679 Interpolation coefficients.
681 Returns
682 -------
683 :class:`numpy.ndarray`
684 Interpolation matrix.
685 """
686 if _is_real(proj):
687 return _backend.idd_reconint(idx + 1, proj)
688 else:
689 return _backend.idz_reconint(idx + 1, proj)
692def reconstruct_skel_matrix(A, k, idx):
693 """
694 Reconstruct skeleton matrix from ID.
696 The skeleton matrix can be reconstructed from the original matrix `A` and its
697 ID rank and indices `k` and `idx`, respectively, as::
699 B = A[:,idx[:k]]
701 The original matrix can then be reconstructed via::
703 numpy.hstack([B, numpy.dot(B, proj)])[:,numpy.argsort(idx)]
705 See also :func:`reconstruct_matrix_from_id` and
706 :func:`reconstruct_interp_matrix`.
708 .. This function automatically detects the matrix data type and calls the
709 appropriate backend. For details, see :func:`_backend.idd_copycols` and
710 :func:`_backend.idz_copycols`.
712 Parameters
713 ----------
714 A : :class:`numpy.ndarray`
715 Original matrix.
716 k : int
717 Rank of ID.
718 idx : :class:`numpy.ndarray`
719 Column index array.
721 Returns
722 -------
723 :class:`numpy.ndarray`
724 Skeleton matrix.
725 """
726 if _is_real(A):
727 return _backend.idd_copycols(A, k, idx + 1)
728 else:
729 return _backend.idz_copycols(A, k, idx + 1)
732def id_to_svd(B, idx, proj):
733 """
734 Convert ID to SVD.
736 The SVD reconstruction of a matrix with skeleton matrix `B` and ID indices and
737 coefficients `idx` and `proj`, respectively, is::
739 U, S, V = id_to_svd(B, idx, proj)
740 A = numpy.dot(U, numpy.dot(numpy.diag(S), V.conj().T))
742 See also :func:`svd`.
744 .. This function automatically detects the matrix data type and calls the
745 appropriate backend. For details, see :func:`_backend.idd_id2svd` and
746 :func:`_backend.idz_id2svd`.
748 Parameters
749 ----------
750 B : :class:`numpy.ndarray`
751 Skeleton matrix.
752 idx : :class:`numpy.ndarray`
753 Column index array.
754 proj : :class:`numpy.ndarray`
755 Interpolation coefficients.
757 Returns
758 -------
759 U : :class:`numpy.ndarray`
760 Left singular vectors.
761 S : :class:`numpy.ndarray`
762 Singular values.
763 V : :class:`numpy.ndarray`
764 Right singular vectors.
765 """
766 if _is_real(B):
767 U, V, S = _backend.idd_id2svd(B, idx + 1, proj)
768 else:
769 U, V, S = _backend.idz_id2svd(B, idx + 1, proj)
770 return U, S, V
773def estimate_spectral_norm(A, its=20):
774 """
775 Estimate spectral norm of a matrix by the randomized power method.
777 .. This function automatically detects the matrix data type and calls the
778 appropriate backend. For details, see :func:`_backend.idd_snorm` and
779 :func:`_backend.idz_snorm`.
781 Parameters
782 ----------
783 A : :class:`scipy.sparse.linalg.LinearOperator`
784 Matrix given as a :class:`scipy.sparse.linalg.LinearOperator` with the
785 `matvec` and `rmatvec` methods (to apply the matrix and its adjoint).
786 its : int, optional
787 Number of power method iterations.
789 Returns
790 -------
791 float
792 Spectral norm estimate.
793 """
794 from scipy.sparse.linalg import aslinearoperator
795 A = aslinearoperator(A)
796 m, n = A.shape
797 def matvec(x):
798 return A.matvec(x)
799 def matveca(x):
800 return A.rmatvec(x)
801 if _is_real(A):
802 return _backend.idd_snorm(m, n, matveca, matvec, its=its)
803 else:
804 return _backend.idz_snorm(m, n, matveca, matvec, its=its)
807def estimate_spectral_norm_diff(A, B, its=20):
808 """
809 Estimate spectral norm of the difference of two matrices by the randomized
810 power method.
812 .. This function automatically detects the matrix data type and calls the
813 appropriate backend. For details, see :func:`_backend.idd_diffsnorm` and
814 :func:`_backend.idz_diffsnorm`.
816 Parameters
817 ----------
818 A : :class:`scipy.sparse.linalg.LinearOperator`
819 First matrix given as a :class:`scipy.sparse.linalg.LinearOperator` with the
820 `matvec` and `rmatvec` methods (to apply the matrix and its adjoint).
821 B : :class:`scipy.sparse.linalg.LinearOperator`
822 Second matrix given as a :class:`scipy.sparse.linalg.LinearOperator` with
823 the `matvec` and `rmatvec` methods (to apply the matrix and its adjoint).
824 its : int, optional
825 Number of power method iterations.
827 Returns
828 -------
829 float
830 Spectral norm estimate of matrix difference.
831 """
832 from scipy.sparse.linalg import aslinearoperator
833 A = aslinearoperator(A)
834 B = aslinearoperator(B)
835 m, n = A.shape
836 def matvec1(x):
837 return A.matvec(x)
838 def matveca1(x):
839 return A.rmatvec(x)
840 def matvec2(x):
841 return B.matvec(x)
842 def matveca2(x):
843 return B.rmatvec(x)
844 if _is_real(A):
845 return _backend.idd_diffsnorm(
846 m, n, matveca1, matveca2, matvec1, matvec2, its=its)
847 else:
848 return _backend.idz_diffsnorm(
849 m, n, matveca1, matveca2, matvec1, matvec2, its=its)
852def svd(A, eps_or_k, rand=True):
853 """
854 Compute SVD of a matrix via an ID.
856 An SVD of a matrix `A` is a factorization::
858 A = numpy.dot(U, numpy.dot(numpy.diag(S), V.conj().T))
860 where `U` and `V` have orthonormal columns and `S` is nonnegative.
862 The SVD can be computed to any relative precision or rank (depending on the
863 value of `eps_or_k`).
865 See also :func:`interp_decomp` and :func:`id_to_svd`.
867 .. This function automatically detects the form of the input parameters and
868 passes them to the appropriate backend. For details, see
869 :func:`_backend.iddp_svd`, :func:`_backend.iddp_asvd`,
870 :func:`_backend.iddp_rsvd`, :func:`_backend.iddr_svd`,
871 :func:`_backend.iddr_asvd`, :func:`_backend.iddr_rsvd`,
872 :func:`_backend.idzp_svd`, :func:`_backend.idzp_asvd`,
873 :func:`_backend.idzp_rsvd`, :func:`_backend.idzr_svd`,
874 :func:`_backend.idzr_asvd`, and :func:`_backend.idzr_rsvd`.
876 Parameters
877 ----------
878 A : :class:`numpy.ndarray` or :class:`scipy.sparse.linalg.LinearOperator`
879 Matrix to be factored, given as either a :class:`numpy.ndarray` or a
880 :class:`scipy.sparse.linalg.LinearOperator` with the `matvec` and
881 `rmatvec` methods (to apply the matrix and its adjoint).
882 eps_or_k : float or int
883 Relative error (if `eps_or_k < 1`) or rank (if `eps_or_k >= 1`) of
884 approximation.
885 rand : bool, optional
886 Whether to use random sampling if `A` is of type :class:`numpy.ndarray`
887 (randomized algorithms are always used if `A` is of type
888 :class:`scipy.sparse.linalg.LinearOperator`).
890 Returns
891 -------
892 U : :class:`numpy.ndarray`
893 Left singular vectors.
894 S : :class:`numpy.ndarray`
895 Singular values.
896 V : :class:`numpy.ndarray`
897 Right singular vectors.
898 """
899 from scipy.sparse.linalg import LinearOperator
901 real = _is_real(A)
903 if isinstance(A, np.ndarray):
904 if eps_or_k < 1:
905 eps = eps_or_k
906 if rand:
907 if real:
908 U, V, S = _backend.iddp_asvd(eps, A)
909 else:
910 if _IS_32BIT:
911 raise _32BIT_ERROR
912 U, V, S = _backend.idzp_asvd(eps, A)
913 else:
914 if real:
915 U, V, S = _backend.iddp_svd(eps, A)
916 else:
917 U, V, S = _backend.idzp_svd(eps, A)
918 else:
919 k = int(eps_or_k)
920 if k > min(A.shape):
921 raise ValueError("Approximation rank {} exceeds min(A.shape) = "
922 " {} ".format(k, min(A.shape)))
923 if rand:
924 if real:
925 U, V, S = _backend.iddr_asvd(A, k)
926 else:
927 if _IS_32BIT:
928 raise _32BIT_ERROR
929 U, V, S = _backend.idzr_asvd(A, k)
930 else:
931 if real:
932 U, V, S = _backend.iddr_svd(A, k)
933 else:
934 U, V, S = _backend.idzr_svd(A, k)
935 elif isinstance(A, LinearOperator):
936 m, n = A.shape
937 def matvec(x):
938 return A.matvec(x)
939 def matveca(x):
940 return A.rmatvec(x)
941 if eps_or_k < 1:
942 eps = eps_or_k
943 if real:
944 U, V, S = _backend.iddp_rsvd(eps, m, n, matveca, matvec)
945 else:
946 if _IS_32BIT:
947 raise _32BIT_ERROR
948 U, V, S = _backend.idzp_rsvd(eps, m, n, matveca, matvec)
949 else:
950 k = int(eps_or_k)
951 if real:
952 U, V, S = _backend.iddr_rsvd(m, n, matveca, matvec, k)
953 else:
954 if _IS_32BIT:
955 raise _32BIT_ERROR
956 U, V, S = _backend.idzr_rsvd(m, n, matveca, matvec, k)
957 else:
958 raise _TYPE_ERROR
959 return U, S, V
962def estimate_rank(A, eps):
963 """
964 Estimate matrix rank to a specified relative precision using randomized
965 methods.
967 The matrix `A` can be given as either a :class:`numpy.ndarray` or a
968 :class:`scipy.sparse.linalg.LinearOperator`, with different algorithms used
969 for each case. If `A` is of type :class:`numpy.ndarray`, then the output
970 rank is typically about 8 higher than the actual numerical rank.
972 .. This function automatically detects the form of the input parameters and
973 passes them to the appropriate backend. For details,
974 see :func:`_backend.idd_estrank`, :func:`_backend.idd_findrank`,
975 :func:`_backend.idz_estrank`, and :func:`_backend.idz_findrank`.
977 Parameters
978 ----------
979 A : :class:`numpy.ndarray` or :class:`scipy.sparse.linalg.LinearOperator`
980 Matrix whose rank is to be estimated, given as either a
981 :class:`numpy.ndarray` or a :class:`scipy.sparse.linalg.LinearOperator`
982 with the `rmatvec` method (to apply the matrix adjoint).
983 eps : float
984 Relative error for numerical rank definition.
986 Returns
987 -------
988 int
989 Estimated matrix rank.
990 """
991 from scipy.sparse.linalg import LinearOperator
993 real = _is_real(A)
995 if isinstance(A, np.ndarray):
996 if real:
997 rank = _backend.idd_estrank(eps, A)
998 else:
999 rank = _backend.idz_estrank(eps, A)
1000 if rank == 0:
1001 # special return value for nearly full rank
1002 rank = min(A.shape)
1003 return rank
1004 elif isinstance(A, LinearOperator):
1005 m, n = A.shape
1006 matveca = A.rmatvec
1007 if real:
1008 return _backend.idd_findrank(eps, m, n, matveca)
1009 else:
1010 return _backend.idz_findrank(eps, m, n, matveca)
1011 else:
1012 raise _TYPE_ERROR