Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/interpolative.py: 12%
181 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 matvec = lambda x: A. matvec(x)
798 matveca = lambda x: A.rmatvec(x)
799 if _is_real(A):
800 return _backend.idd_snorm(m, n, matveca, matvec, its=its)
801 else:
802 return _backend.idz_snorm(m, n, matveca, matvec, its=its)
805def estimate_spectral_norm_diff(A, B, its=20):
806 """
807 Estimate spectral norm of the difference of two matrices by the randomized
808 power method.
810 .. This function automatically detects the matrix data type and calls the
811 appropriate backend. For details, see :func:`_backend.idd_diffsnorm` and
812 :func:`_backend.idz_diffsnorm`.
814 Parameters
815 ----------
816 A : :class:`scipy.sparse.linalg.LinearOperator`
817 First matrix given as a :class:`scipy.sparse.linalg.LinearOperator` with the
818 `matvec` and `rmatvec` methods (to apply the matrix and its adjoint).
819 B : :class:`scipy.sparse.linalg.LinearOperator`
820 Second matrix given as a :class:`scipy.sparse.linalg.LinearOperator` with
821 the `matvec` and `rmatvec` methods (to apply the matrix and its adjoint).
822 its : int, optional
823 Number of power method iterations.
825 Returns
826 -------
827 float
828 Spectral norm estimate of matrix difference.
829 """
830 from scipy.sparse.linalg import aslinearoperator
831 A = aslinearoperator(A)
832 B = aslinearoperator(B)
833 m, n = A.shape
834 matvec1 = lambda x: A. matvec(x)
835 matveca1 = lambda x: A.rmatvec(x)
836 matvec2 = lambda x: B. matvec(x)
837 matveca2 = lambda x: B.rmatvec(x)
838 if _is_real(A):
839 return _backend.idd_diffsnorm(
840 m, n, matveca1, matveca2, matvec1, matvec2, its=its)
841 else:
842 return _backend.idz_diffsnorm(
843 m, n, matveca1, matveca2, matvec1, matvec2, its=its)
846def svd(A, eps_or_k, rand=True):
847 """
848 Compute SVD of a matrix via an ID.
850 An SVD of a matrix `A` is a factorization::
852 A = numpy.dot(U, numpy.dot(numpy.diag(S), V.conj().T))
854 where `U` and `V` have orthonormal columns and `S` is nonnegative.
856 The SVD can be computed to any relative precision or rank (depending on the
857 value of `eps_or_k`).
859 See also :func:`interp_decomp` and :func:`id_to_svd`.
861 .. This function automatically detects the form of the input parameters and
862 passes them to the appropriate backend. For details, see
863 :func:`_backend.iddp_svd`, :func:`_backend.iddp_asvd`,
864 :func:`_backend.iddp_rsvd`, :func:`_backend.iddr_svd`,
865 :func:`_backend.iddr_asvd`, :func:`_backend.iddr_rsvd`,
866 :func:`_backend.idzp_svd`, :func:`_backend.idzp_asvd`,
867 :func:`_backend.idzp_rsvd`, :func:`_backend.idzr_svd`,
868 :func:`_backend.idzr_asvd`, and :func:`_backend.idzr_rsvd`.
870 Parameters
871 ----------
872 A : :class:`numpy.ndarray` or :class:`scipy.sparse.linalg.LinearOperator`
873 Matrix to be factored, given as either a :class:`numpy.ndarray` or a
874 :class:`scipy.sparse.linalg.LinearOperator` with the `matvec` and
875 `rmatvec` methods (to apply the matrix and its adjoint).
876 eps_or_k : float or int
877 Relative error (if `eps_or_k < 1`) or rank (if `eps_or_k >= 1`) of
878 approximation.
879 rand : bool, optional
880 Whether to use random sampling if `A` is of type :class:`numpy.ndarray`
881 (randomized algorithms are always used if `A` is of type
882 :class:`scipy.sparse.linalg.LinearOperator`).
884 Returns
885 -------
886 U : :class:`numpy.ndarray`
887 Left singular vectors.
888 S : :class:`numpy.ndarray`
889 Singular values.
890 V : :class:`numpy.ndarray`
891 Right singular vectors.
892 """
893 from scipy.sparse.linalg import LinearOperator
895 real = _is_real(A)
897 if isinstance(A, np.ndarray):
898 if eps_or_k < 1:
899 eps = eps_or_k
900 if rand:
901 if real:
902 U, V, S = _backend.iddp_asvd(eps, A)
903 else:
904 if _IS_32BIT:
905 raise _32BIT_ERROR
906 U, V, S = _backend.idzp_asvd(eps, A)
907 else:
908 if real:
909 U, V, S = _backend.iddp_svd(eps, A)
910 else:
911 U, V, S = _backend.idzp_svd(eps, A)
912 else:
913 k = int(eps_or_k)
914 if k > min(A.shape):
915 raise ValueError("Approximation rank %s exceeds min(A.shape) = "
916 " %s " % (k, min(A.shape)))
917 if rand:
918 if real:
919 U, V, S = _backend.iddr_asvd(A, k)
920 else:
921 if _IS_32BIT:
922 raise _32BIT_ERROR
923 U, V, S = _backend.idzr_asvd(A, k)
924 else:
925 if real:
926 U, V, S = _backend.iddr_svd(A, k)
927 else:
928 U, V, S = _backend.idzr_svd(A, k)
929 elif isinstance(A, LinearOperator):
930 m, n = A.shape
931 matvec = lambda x: A.matvec(x)
932 matveca = lambda x: A.rmatvec(x)
933 if eps_or_k < 1:
934 eps = eps_or_k
935 if real:
936 U, V, S = _backend.iddp_rsvd(eps, m, n, matveca, matvec)
937 else:
938 if _IS_32BIT:
939 raise _32BIT_ERROR
940 U, V, S = _backend.idzp_rsvd(eps, m, n, matveca, matvec)
941 else:
942 k = int(eps_or_k)
943 if real:
944 U, V, S = _backend.iddr_rsvd(m, n, matveca, matvec, k)
945 else:
946 if _IS_32BIT:
947 raise _32BIT_ERROR
948 U, V, S = _backend.idzr_rsvd(m, n, matveca, matvec, k)
949 else:
950 raise _TYPE_ERROR
951 return U, S, V
954def estimate_rank(A, eps):
955 """
956 Estimate matrix rank to a specified relative precision using randomized
957 methods.
959 The matrix `A` can be given as either a :class:`numpy.ndarray` or a
960 :class:`scipy.sparse.linalg.LinearOperator`, with different algorithms used
961 for each case. If `A` is of type :class:`numpy.ndarray`, then the output
962 rank is typically about 8 higher than the actual numerical rank.
964 .. This function automatically detects the form of the input parameters and
965 passes them to the appropriate backend. For details,
966 see :func:`_backend.idd_estrank`, :func:`_backend.idd_findrank`,
967 :func:`_backend.idz_estrank`, and :func:`_backend.idz_findrank`.
969 Parameters
970 ----------
971 A : :class:`numpy.ndarray` or :class:`scipy.sparse.linalg.LinearOperator`
972 Matrix whose rank is to be estimated, given as either a
973 :class:`numpy.ndarray` or a :class:`scipy.sparse.linalg.LinearOperator`
974 with the `rmatvec` method (to apply the matrix adjoint).
975 eps : float
976 Relative error for numerical rank definition.
978 Returns
979 -------
980 int
981 Estimated matrix rank.
982 """
983 from scipy.sparse.linalg import LinearOperator
985 real = _is_real(A)
987 if isinstance(A, np.ndarray):
988 if real:
989 rank = _backend.idd_estrank(eps, A)
990 else:
991 rank = _backend.idz_estrank(eps, A)
992 if rank == 0:
993 # special return value for nearly full rank
994 rank = min(A.shape)
995 return rank
996 elif isinstance(A, LinearOperator):
997 m, n = A.shape
998 matveca = A.rmatvec
999 if real:
1000 return _backend.idd_findrank(eps, m, n, matveca)
1001 else:
1002 return _backend.idz_findrank(eps, m, n, matveca)
1003 else:
1004 raise _TYPE_ERROR