Why Gemfury?
Push, build, and install
RubyGems
npm packages
Python packages
Maven artifacts
PHP packages
Go Modules
Bower components
Debian packages
RPM packages
NuGet packages
/ linalg / _interpolative_backend.py
#******************************************************************************
# Copyright (C) 2013 Kenneth L. Ho
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimer. Redistributions in binary
# form must reproduce the above copyright notice, this list of conditions and
# the following disclaimer in the documentation and/or other materials
# provided with the distribution.
#
# None of the names of the copyright holders may be used to endorse or
# promote products derived from this software without specific prior written
# permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#******************************************************************************
"""
Direct wrappers for Fortran `id_dist` backend.
"""
import scipy.linalg._interpolative as _id
import numpy as np
_RETCODE_ERROR = RuntimeError("nonzero return code")
#------------------------------------------------------------------------------
# id_rand.f
#------------------------------------------------------------------------------
def id_srand(n):
"""
Generate standard uniform pseudorandom numbers via a very efficient lagged
Fibonacci method.
:param n:
Number of pseudorandom numbers to generate.
:type n: int
:return:
Pseudorandom numbers.
:rtype: :class:`numpy.ndarray`
"""
return _id.id_srand(n)
def id_srandi(t):
"""
Initialize seed values for :func:`id_srand` (any appropriately random
numbers will do).
:param t:
Array of 55 seed values.
:type t: :class:`numpy.ndarray`
"""
t = np.asfortranarray(t)
_id.id_srandi(t)
def id_srando():
"""
Reset seed values to their original values.
"""
_id.id_srando()
#------------------------------------------------------------------------------
# idd_frm.f
#------------------------------------------------------------------------------
def idd_frm(n, w, x):
"""
Transform real vector via a composition of Rokhlin's random transform,
random subselection, and an FFT.
In contrast to :func:`idd_sfrm`, this routine works best when the length of
the transformed vector is the power-of-two integer output by
:func:`idd_frmi`, or when the length is not specified but instead
determined a posteriori from the output. The returned transformed vector is
randomly permuted.
:param n:
Greatest power-of-two integer satisfying `n <= x.size` as obtained from
:func:`idd_frmi`; `n` is also the length of the output vector.
:type n: int
:param w:
Initialization array constructed by :func:`idd_frmi`.
:type w: :class:`numpy.ndarray`
:param x:
Vector to be transformed.
:type x: :class:`numpy.ndarray`
:return:
Transformed vector.
:rtype: :class:`numpy.ndarray`
"""
return _id.idd_frm(n, w, x)
def idd_sfrm(l, n, w, x):
"""
Transform real vector via a composition of Rokhlin's random transform,
random subselection, and an FFT.
In contrast to :func:`idd_frm`, this routine works best when the length of
the transformed vector is known a priori.
:param l:
Length of transformed vector, satisfying `l <= n`.
:type l: int
:param n:
Greatest power-of-two integer satisfying `n <= x.size` as obtained from
:func:`idd_sfrmi`.
:type n: int
:param w:
Initialization array constructed by :func:`idd_sfrmi`.
:type w: :class:`numpy.ndarray`
:param x:
Vector to be transformed.
:type x: :class:`numpy.ndarray`
:return:
Transformed vector.
:rtype: :class:`numpy.ndarray`
"""
return _id.idd_sfrm(l, n, w, x)
def idd_frmi(m):
"""
Initialize data for :func:`idd_frm`.
:param m:
Length of vector to be transformed.
:type m: int
:return:
Greatest power-of-two integer `n` satisfying `n <= m`.
:rtype: int
:return:
Initialization array to be used by :func:`idd_frm`.
:rtype: :class:`numpy.ndarray`
"""
return _id.idd_frmi(m)
def idd_sfrmi(l, m):
"""
Initialize data for :func:`idd_sfrm`.
:param l:
Length of output transformed vector.
:type l: int
:param m:
Length of the vector to be transformed.
:type m: int
:return:
Greatest power-of-two integer `n` satisfying `n <= m`.
:rtype: int
:return:
Initialization array to be used by :func:`idd_sfrm`.
:rtype: :class:`numpy.ndarray`
"""
return _id.idd_sfrmi(l, m)
#------------------------------------------------------------------------------
# idd_id.f
#------------------------------------------------------------------------------
def iddp_id(eps, A):
"""
Compute ID of a real matrix to a specified relative precision.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
k, idx, rnorms = _id.iddp_id(eps, A)
n = A.shape[1]
proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def iddr_id(A, k):
"""
Compute ID of a real matrix to a specified rank.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
idx, rnorms = _id.iddr_id(A, k)
n = A.shape[1]
proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
return idx, proj
def idd_reconid(B, idx, proj):
"""
Reconstruct matrix from real ID.
:param B:
Skeleton matrix.
:type B: :class:`numpy.ndarray`
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Reconstructed matrix.
:rtype: :class:`numpy.ndarray`
"""
B = np.asfortranarray(B)
if proj.size > 0:
return _id.idd_reconid(B, idx, proj)
else:
return B[:, np.argsort(idx)]
def idd_reconint(idx, proj):
"""
Reconstruct interpolation matrix from real ID.
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Interpolation matrix.
:rtype: :class:`numpy.ndarray`
"""
return _id.idd_reconint(idx, proj)
def idd_copycols(A, k, idx):
"""
Reconstruct skeleton matrix from real ID.
:param A:
Original matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:return:
Skeleton matrix.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
return _id.idd_copycols(A, k, idx)
#------------------------------------------------------------------------------
# idd_id2svd.f
#------------------------------------------------------------------------------
def idd_id2svd(B, idx, proj):
"""
Convert real ID to SVD.
:param B:
Skeleton matrix.
:type B: :class:`numpy.ndarray`
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
B = np.asfortranarray(B)
U, V, S, ier = _id.idd_id2svd(B, idx, proj)
if ier:
raise _RETCODE_ERROR
return U, V, S
#------------------------------------------------------------------------------
# idd_snorm.f
#------------------------------------------------------------------------------
def idd_snorm(m, n, matvect, matvec, its=20):
"""
Estimate spectral norm of a real matrix by the randomized power method.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param its:
Number of power method iterations.
:type its: int
:return:
Spectral norm estimate.
:rtype: float
"""
snorm, v = _id.idd_snorm(m, n, matvect, matvec, its)
return snorm
def idd_diffsnorm(m, n, matvect, matvect2, matvec, matvec2, its=20):
"""
Estimate spectral norm of the difference of two real matrices by the
randomized power method.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the transpose of the first matrix to a vector, with
call signature `y = matvect(x)`, where `x` and `y` are the input and
output vectors, respectively.
:type matvect: function
:param matvect2:
Function to apply the transpose of the second matrix to a vector, with
call signature `y = matvect2(x)`, where `x` and `y` are the input and
output vectors, respectively.
:type matvect2: function
:param matvec:
Function to apply the first matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param matvec2:
Function to apply the second matrix to a vector, with call signature
`y = matvec2(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec2: function
:param its:
Number of power method iterations.
:type its: int
:return:
Spectral norm estimate of matrix difference.
:rtype: float
"""
return _id.idd_diffsnorm(m, n, matvect, matvect2, matvec, matvec2, its)
#------------------------------------------------------------------------------
# idd_svd.f
#------------------------------------------------------------------------------
def iddr_svd(A, k):
"""
Compute SVD of a real matrix to a specified rank.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
U, V, S, ier = _id.iddr_svd(A, k)
if ier:
raise _RETCODE_ERROR
return U, V, S
def iddp_svd(eps, A):
"""
Compute SVD of a real matrix to a specified relative precision.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
k, iU, iV, iS, w, ier = _id.iddp_svd(eps, A)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# iddp_aid.f
#------------------------------------------------------------------------------
def iddp_aid(eps, A):
"""
Compute ID of a real matrix to a specified relative precision using random
sampling.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, w = idd_frmi(m)
proj = np.empty(n*(2*n2 + 1) + n2 + 1, order='F')
k, idx, proj = _id.iddp_aid(eps, A, w, proj)
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def idd_estrank(eps, A):
"""
Estimate rank of a real matrix to a specified relative precision using
random sampling.
The output rank is typically about 8 higher than the actual rank.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank estimate.
:rtype: int
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, w = idd_frmi(m)
ra = np.empty(n*n2 + (n + 1)*(n2 + 1), order='F')
k, ra = _id.idd_estrank(eps, A, w, ra)
return k
#------------------------------------------------------------------------------
# iddp_asvd.f
#------------------------------------------------------------------------------
def iddp_asvd(eps, A):
"""
Compute SVD of a real matrix to a specified relative precision using random
sampling.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, winit = _id.idd_frmi(m)
w = np.empty(
max((min(m, n) + 1)*(3*m + 5*n + 1) + 25*min(m, n)**2,
(2*n + 1)*(n2 + 1)),
order='F')
k, iU, iV, iS, w, ier = _id.iddp_asvd(eps, A, winit, w)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# iddp_rid.f
#------------------------------------------------------------------------------
def iddp_rid(eps, m, n, matvect):
"""
Compute ID of a real matrix to a specified relative precision using random
matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
proj = np.empty(m + 1 + 2*n*(min(m, n) + 1), order='F')
k, idx, proj, ier = _id.iddp_rid(eps, m, n, matvect, proj)
if ier != 0:
raise _RETCODE_ERROR
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def idd_findrank(eps, m, n, matvect):
"""
Estimate rank of a real matrix to a specified relative precision using
random matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:return:
Rank estimate.
:rtype: int
"""
k, ra, ier = _id.idd_findrank(eps, m, n, matvect)
if ier:
raise _RETCODE_ERROR
return k
#------------------------------------------------------------------------------
# iddp_rsvd.f
#------------------------------------------------------------------------------
def iddp_rsvd(eps, m, n, matvect, matvec):
"""
Compute SVD of a real matrix to a specified relative precision using random
matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
k, iU, iV, iS, w, ier = _id.iddp_rsvd(eps, m, n, matvect, matvec)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# iddr_aid.f
#------------------------------------------------------------------------------
def iddr_aid(A, k):
"""
Compute ID of a real matrix to a specified rank using random sampling.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
w = iddr_aidi(m, n, k)
idx, proj = _id.iddr_aid(A, k, w)
if k == n:
proj = np.empty((k, n-k), dtype='float64', order='F')
else:
proj = proj.reshape((k, n-k), order='F')
return idx, proj
def iddr_aidi(m, n, k):
"""
Initialize array for :func:`iddr_aid`.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param k:
Rank of ID.
:type k: int
:return:
Initialization array to be used by :func:`iddr_aid`.
:rtype: :class:`numpy.ndarray`
"""
return _id.iddr_aidi(m, n, k)
#------------------------------------------------------------------------------
# iddr_asvd.f
#------------------------------------------------------------------------------
def iddr_asvd(A, k):
"""
Compute SVD of a real matrix to a specified rank using random sampling.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
w = np.empty((2*k + 28)*m + (6*k + 21)*n + 25*k**2 + 100, order='F')
w_ = iddr_aidi(m, n, k)
w[:w_.size] = w_
U, V, S, ier = _id.iddr_asvd(A, k, w)
if ier != 0:
raise _RETCODE_ERROR
return U, V, S
#------------------------------------------------------------------------------
# iddr_rid.f
#------------------------------------------------------------------------------
def iddr_rid(m, n, matvect, k):
"""
Compute ID of a real matrix to a specified rank using random matrix-vector
multiplication.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
idx, proj = _id.iddr_rid(m, n, matvect, k)
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return idx, proj
#------------------------------------------------------------------------------
# iddr_rsvd.f
#------------------------------------------------------------------------------
def iddr_rsvd(m, n, matvect, matvec, k):
"""
Compute SVD of a real matrix to a specified rank using random matrix-vector
multiplication.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matvect:
Function to apply the matrix transpose to a vector, with call signature
`y = matvect(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvect: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
U, V, S, ier = _id.iddr_rsvd(m, n, matvect, matvec, k)
if ier != 0:
raise _RETCODE_ERROR
return U, V, S
#------------------------------------------------------------------------------
# idz_frm.f
#------------------------------------------------------------------------------
def idz_frm(n, w, x):
"""
Transform complex vector via a composition of Rokhlin's random transform,
random subselection, and an FFT.
In contrast to :func:`idz_sfrm`, this routine works best when the length of
the transformed vector is the power-of-two integer output by
:func:`idz_frmi`, or when the length is not specified but instead
determined a posteriori from the output. The returned transformed vector is
randomly permuted.
:param n:
Greatest power-of-two integer satisfying `n <= x.size` as obtained from
:func:`idz_frmi`; `n` is also the length of the output vector.
:type n: int
:param w:
Initialization array constructed by :func:`idz_frmi`.
:type w: :class:`numpy.ndarray`
:param x:
Vector to be transformed.
:type x: :class:`numpy.ndarray`
:return:
Transformed vector.
:rtype: :class:`numpy.ndarray`
"""
return _id.idz_frm(n, w, x)
def idz_sfrm(l, n, w, x):
"""
Transform complex vector via a composition of Rokhlin's random transform,
random subselection, and an FFT.
In contrast to :func:`idz_frm`, this routine works best when the length of
the transformed vector is known a priori.
:param l:
Length of transformed vector, satisfying `l <= n`.
:type l: int
:param n:
Greatest power-of-two integer satisfying `n <= x.size` as obtained from
:func:`idz_sfrmi`.
:type n: int
:param w:
Initialization array constructed by :func:`idd_sfrmi`.
:type w: :class:`numpy.ndarray`
:param x:
Vector to be transformed.
:type x: :class:`numpy.ndarray`
:return:
Transformed vector.
:rtype: :class:`numpy.ndarray`
"""
return _id.idz_sfrm(l, n, w, x)
def idz_frmi(m):
"""
Initialize data for :func:`idz_frm`.
:param m:
Length of vector to be transformed.
:type m: int
:return:
Greatest power-of-two integer `n` satisfying `n <= m`.
:rtype: int
:return:
Initialization array to be used by :func:`idz_frm`.
:rtype: :class:`numpy.ndarray`
"""
return _id.idz_frmi(m)
def idz_sfrmi(l, m):
"""
Initialize data for :func:`idz_sfrm`.
:param l:
Length of output transformed vector.
:type l: int
:param m:
Length of the vector to be transformed.
:type m: int
:return:
Greatest power-of-two integer `n` satisfying `n <= m`.
:rtype: int
:return:
Initialization array to be used by :func:`idz_sfrm`.
:rtype: :class:`numpy.ndarray`
"""
return _id.idz_sfrmi(l, m)
#------------------------------------------------------------------------------
# idz_id.f
#------------------------------------------------------------------------------
def idzp_id(eps, A):
"""
Compute ID of a complex matrix to a specified relative precision.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
k, idx, rnorms = _id.idzp_id(eps, A)
n = A.shape[1]
proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def idzr_id(A, k):
"""
Compute ID of a complex matrix to a specified rank.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
idx, rnorms = _id.idzr_id(A, k)
n = A.shape[1]
proj = A.T.ravel()[:k*(n-k)].reshape((k, n-k), order='F')
return idx, proj
def idz_reconid(B, idx, proj):
"""
Reconstruct matrix from complex ID.
:param B:
Skeleton matrix.
:type B: :class:`numpy.ndarray`
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Reconstructed matrix.
:rtype: :class:`numpy.ndarray`
"""
B = np.asfortranarray(B)
if proj.size > 0:
return _id.idz_reconid(B, idx, proj)
else:
return B[:, np.argsort(idx)]
def idz_reconint(idx, proj):
"""
Reconstruct interpolation matrix from complex ID.
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Interpolation matrix.
:rtype: :class:`numpy.ndarray`
"""
return _id.idz_reconint(idx, proj)
def idz_copycols(A, k, idx):
"""
Reconstruct skeleton matrix from complex ID.
:param A:
Original matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:return:
Skeleton matrix.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
return _id.idz_copycols(A, k, idx)
#------------------------------------------------------------------------------
# idz_id2svd.f
#------------------------------------------------------------------------------
def idz_id2svd(B, idx, proj):
"""
Convert complex ID to SVD.
:param B:
Skeleton matrix.
:type B: :class:`numpy.ndarray`
:param idx:
Column index array.
:type idx: :class:`numpy.ndarray`
:param proj:
Interpolation coefficients.
:type proj: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
B = np.asfortranarray(B)
U, V, S, ier = _id.idz_id2svd(B, idx, proj)
if ier:
raise _RETCODE_ERROR
return U, V, S
#------------------------------------------------------------------------------
# idz_snorm.f
#------------------------------------------------------------------------------
def idz_snorm(m, n, matveca, matvec, its=20):
"""
Estimate spectral norm of a complex matrix by the randomized power method.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param its:
Number of power method iterations.
:type its: int
:return:
Spectral norm estimate.
:rtype: float
"""
snorm, v = _id.idz_snorm(m, n, matveca, matvec, its)
return snorm
def idz_diffsnorm(m, n, matveca, matveca2, matvec, matvec2, its=20):
"""
Estimate spectral norm of the difference of two complex matrices by the
randomized power method.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the adjoint of the first matrix to a vector, with
call signature `y = matveca(x)`, where `x` and `y` are the input and
output vectors, respectively.
:type matveca: function
:param matveca2:
Function to apply the adjoint of the second matrix to a vector, with
call signature `y = matveca2(x)`, where `x` and `y` are the input and
output vectors, respectively.
:type matveca2: function
:param matvec:
Function to apply the first matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param matvec2:
Function to apply the second matrix to a vector, with call signature
`y = matvec2(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec2: function
:param its:
Number of power method iterations.
:type its: int
:return:
Spectral norm estimate of matrix difference.
:rtype: float
"""
return _id.idz_diffsnorm(m, n, matveca, matveca2, matvec, matvec2, its)
#------------------------------------------------------------------------------
# idz_svd.f
#------------------------------------------------------------------------------
def idzr_svd(A, k):
"""
Compute SVD of a complex matrix to a specified rank.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
U, V, S, ier = _id.idzr_svd(A, k)
if ier:
raise _RETCODE_ERROR
return U, V, S
def idzp_svd(eps, A):
"""
Compute SVD of a complex matrix to a specified relative precision.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
k, iU, iV, iS, w, ier = _id.idzp_svd(eps, A)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# idzp_aid.f
#------------------------------------------------------------------------------
def idzp_aid(eps, A):
"""
Compute ID of a complex matrix to a specified relative precision using
random sampling.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, w = idz_frmi(m)
proj = np.empty(n*(2*n2 + 1) + n2 + 1, dtype='complex128', order='F')
k, idx, proj = _id.idzp_aid(eps, A, w, proj)
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def idz_estrank(eps, A):
"""
Estimate rank of a complex matrix to a specified relative precision using
random sampling.
The output rank is typically about 8 higher than the actual rank.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Rank estimate.
:rtype: int
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, w = idz_frmi(m)
ra = np.empty(n*n2 + (n + 1)*(n2 + 1), dtype='complex128', order='F')
k, ra = _id.idz_estrank(eps, A, w, ra)
return k
#------------------------------------------------------------------------------
# idzp_asvd.f
#------------------------------------------------------------------------------
def idzp_asvd(eps, A):
"""
Compute SVD of a complex matrix to a specified relative precision using
random sampling.
:param eps:
Relative precision.
:type eps: float
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
n2, winit = _id.idz_frmi(m)
w = np.empty(
max((min(m, n) + 1)*(3*m + 5*n + 11) + 8*min(m, n)**2,
(2*n + 1)*(n2 + 1)),
dtype=np.complex128, order='F')
k, iU, iV, iS, w, ier = _id.idzp_asvd(eps, A, winit, w)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# idzp_rid.f
#------------------------------------------------------------------------------
def idzp_rid(eps, m, n, matveca):
"""
Compute ID of a complex matrix to a specified relative precision using
random matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:return:
Rank of ID.
:rtype: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
proj = np.empty(
m + 1 + 2*n*(min(m, n) + 1),
dtype=np.complex128, order='F')
k, idx, proj, ier = _id.idzp_rid(eps, m, n, matveca, proj)
if ier:
raise _RETCODE_ERROR
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return k, idx, proj
def idz_findrank(eps, m, n, matveca):
"""
Estimate rank of a complex matrix to a specified relative precision using
random matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:return:
Rank estimate.
:rtype: int
"""
k, ra, ier = _id.idz_findrank(eps, m, n, matveca)
if ier:
raise _RETCODE_ERROR
return k
#------------------------------------------------------------------------------
# idzp_rsvd.f
#------------------------------------------------------------------------------
def idzp_rsvd(eps, m, n, matveca, matvec):
"""
Compute SVD of a complex matrix to a specified relative precision using
random matrix-vector multiplication.
:param eps:
Relative precision.
:type eps: float
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
k, iU, iV, iS, w, ier = _id.idzp_rsvd(eps, m, n, matveca, matvec)
if ier:
raise _RETCODE_ERROR
U = w[iU-1:iU+m*k-1].reshape((m, k), order='F')
V = w[iV-1:iV+n*k-1].reshape((n, k), order='F')
S = w[iS-1:iS+k-1]
return U, V, S
#------------------------------------------------------------------------------
# idzr_aid.f
#------------------------------------------------------------------------------
def idzr_aid(A, k):
"""
Compute ID of a complex matrix to a specified rank using random sampling.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
w = idzr_aidi(m, n, k)
idx, proj = _id.idzr_aid(A, k, w)
if k == n:
proj = np.empty((k, n-k), dtype='complex128', order='F')
else:
proj = proj.reshape((k, n-k), order='F')
return idx, proj
def idzr_aidi(m, n, k):
"""
Initialize array for :func:`idzr_aid`.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param k:
Rank of ID.
:type k: int
:return:
Initialization array to be used by :func:`idzr_aid`.
:rtype: :class:`numpy.ndarray`
"""
return _id.idzr_aidi(m, n, k)
#------------------------------------------------------------------------------
# idzr_asvd.f
#------------------------------------------------------------------------------
def idzr_asvd(A, k):
"""
Compute SVD of a complex matrix to a specified rank using random sampling.
:param A:
Matrix.
:type A: :class:`numpy.ndarray`
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
A = np.asfortranarray(A)
m, n = A.shape
w = np.empty(
(2*k + 22)*m + (6*k + 21)*n + 8*k**2 + 10*k + 90,
dtype='complex128', order='F')
w_ = idzr_aidi(m, n, k)
w[:w_.size] = w_
U, V, S, ier = _id.idzr_asvd(A, k, w)
if ier:
raise _RETCODE_ERROR
return U, V, S
#------------------------------------------------------------------------------
# idzr_rid.f
#------------------------------------------------------------------------------
def idzr_rid(m, n, matveca, k):
"""
Compute ID of a complex matrix to a specified rank using random
matrix-vector multiplication.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:param k:
Rank of ID.
:type k: int
:return:
Column index array.
:rtype: :class:`numpy.ndarray`
:return:
Interpolation coefficients.
:rtype: :class:`numpy.ndarray`
"""
idx, proj = _id.idzr_rid(m, n, matveca, k)
proj = proj[:k*(n-k)].reshape((k, n-k), order='F')
return idx, proj
#------------------------------------------------------------------------------
# idzr_rsvd.f
#------------------------------------------------------------------------------
def idzr_rsvd(m, n, matveca, matvec, k):
"""
Compute SVD of a complex matrix to a specified rank using random
matrix-vector multiplication.
:param m:
Matrix row dimension.
:type m: int
:param n:
Matrix column dimension.
:type n: int
:param matveca:
Function to apply the matrix adjoint to a vector, with call signature
`y = matveca(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matveca: function
:param matvec:
Function to apply the matrix to a vector, with call signature
`y = matvec(x)`, where `x` and `y` are the input and output vectors,
respectively.
:type matvec: function
:param k:
Rank of SVD.
:type k: int
:return:
Left singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Right singular vectors.
:rtype: :class:`numpy.ndarray`
:return:
Singular values.
:rtype: :class:`numpy.ndarray`
"""
U, V, S, ier = _id.idzr_rsvd(m, n, matveca, matvec, k)
if ier:
raise _RETCODE_ERROR
return U, V, S