Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixes #7 via sketch module. #19

Merged
merged 7 commits into from
May 26, 2018
Merged

Fixes #7 via sketch module. #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
377e997
WIP: working on sketch module\nWorking on writing a seperate sketch m…
jknox13 May 18, 2018
77c8cef
WIP: need to reconcile edits with package.\nWork from the plane, need…
jknox13 May 18, 2018
1838bb8
WIP: Tests passing for sketch module.
jknox13 May 21, 2018
89b7030
WIP: Finished phasing out old sketch module.
jknox13 May 23, 2018
0b4815a
WIP: change parameter names, need to fix tests.
jknox13 May 23, 2018
e4deccd
Tests now align with PEP parameter changes.
jknox13 May 24, 2018
c685abd
Move calling of subspace iterations.
jknox13 May 24, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 28 additions & 22 deletions ristretto/cur.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
from .interp_decomp import interp_decomp, rinterp_decomp


def cur(A, k=None, index_set=False):
def cur(A, rank=None, index_set=False):
"""CUR decomposition.

Algorithm for computing the low-rank CUR
Expand All @@ -28,10 +28,10 @@ def cur(A, k=None, index_set=False):
Parameters
----------
A : array_like, shape `(m, n)`.
Input matrix.
Input array.

k : integer, `k << min{m,n}`.
Target rank.
rank : integer
Target rank. Best if `rank << min{m,n}`

index_set: str `{'True', 'False'}`, default: `index_set='False'`.
'True' : Return column/row index set instead of `C` and `R`.
Expand All @@ -57,13 +57,13 @@ def cur(A, k=None, index_set=False):
(available at `arXiv <http://arxiv.org/abs/1502.05366>`_).
"""
# compute column ID
J, V = interp_decomp(A, k=k, mode='column', index_set=True)
J, V = interp_decomp(A, rank, mode='column', index_set=True)

# select column subset
C = A[:, J]

# compute row ID of C
Z, I = interp_decomp(C, k=k, mode='row', index_set=True)
Z, I = interp_decomp(C, rank, mode='row', index_set=True)

# select row subset
R = A[I, :]
Expand All @@ -77,31 +77,37 @@ def cur(A, k=None, index_set=False):
return C, U, R


def rcur(A, k=None, p=10, q=1, index_set=False, random_state=None):
def rcur(A, rank, oversample=10, n_subspace=2, index_set=False, random_state=None):
"""Randomized CUR decomposition.

Randomized algorithm for computing the approximate low-rank CUR
decomposition of a rectangular `(m, n)` matrix `A`, with target rank `k << min{m, n}`.
decomposition of a rectangular `(m, n)` matrix `A`, with target rank `rank << min{m, n}`.
Input matrix is factored as `A = C * U * R`, using the column/row pivoted QR decomposition.
The factor matrix `C` is formed of a subset of columns of `A`,
also called the partial column skeleton. The factor matrix `R` is formed as
a subset of rows of `A` also called the partial row skeleton.
The factor matrix `U` is formed so that `U = C**-1 * A * R**-1` is satisfied.

The quality of the approximation can be controlled via the oversampling
parameter `oversample` and `n_subspace` which specifies the number of
subspace iterations.


Parameters
----------
A : array_like, shape `(m, n)`.
Input matrix.
Input array.

k : integer, `k << min{m,n}`.
Target rank.
rank : integer
Target rank. Best if `rank << min{m,n}`

p : integer, default: `p=10`.
Parameter to control oversampling.
oversample : integer, optional (default: 10)
Controls the oversampling of column space. Increasing this parameter
may improve numerical accuracy.

q : integer, default: `q=1`.
Parameter to control number of power (subspace) iterations.
n_subspace : integer, default: 2.
Parameter to control number of subspace iterations. Increasing this
parameter may improve numerical accuracy.

index_set: str `{'True', 'False'}`, default: `index_set='False'`.
'True' : Return column/row index set instead of `C` and `R`.
Expand All @@ -114,13 +120,13 @@ def rcur(A, k=None, p=10, q=1, index_set=False, random_state=None):

Returns
-------
C: array_like, shape `(m, k)`.
C: array_like, shape `(m, rank)`.
Partial column skeleton.

U : array_like, shape `(k, k)`.
U : array_like, shape `(rank, rank)`.
Well-conditioned matrix.

R : array_like, shape `(k, n)`.
R : array_like, shape `(rank, n)`.
Partial row skeleton.


Expand All @@ -133,15 +139,15 @@ def rcur(A, k=None, p=10, q=1, index_set=False, random_state=None):
(available at `arXiv <http://arxiv.org/abs/1502.05366>`_).
"""
# Compute column ID
J, V = rinterp_decomp(A, k=k, p=p, q=q, mode='column', index_set=True,
random_state=random_state)
J, V = rinterp_decomp(A, rank, oversample=oversample, n_subspace=n_subspace,
mode='column', index_set=True, random_state=random_state)

# Select column subset
C = A[:, J]

# Compute row ID of C
Z, I = rinterp_decomp(C, k=k, p=p, q=q, mode='row', index_set=True,
random_state=random_state)
Z, I = rinterp_decomp(A, rank, oversample=oversample, n_subspace=n_subspace,
mode='row', index_set=True, random_state=random_state)

# Select row subset
R = A[I, :]
Expand Down
81 changes: 37 additions & 44 deletions ristretto/dmd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ def _get_amplitudes(F, A):
return result[0]


def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,
def dmd(A, rank=None, dt=1, modes='exact', return_amplitudes=False,
return_vandermonde=False, order=True):
"""Dynamic Mode Decomposition.

Expand All @@ -36,14 +36,14 @@ def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,

Parameters
----------
A : array_like
Real/complex input matrix `a` with dimensions `(m, n)`.
A : array_like, shape `(m, n)`.
Input array.

dt : scalar or array_like
Factor specifying the time difference between the observations.
rank : int
If `rank < (n-1)` low-rank Dynamic Mode Decomposition is computed.

k : int, optional
If `k < (n-1)` low-rank Dynamic Mode Decomposition is computed.
dt : scalar or array_like, optional (default: 1)
Factor specifying the time difference between the observations.

modes : str `{'standard', 'exact', 'exact_scaled'}`
- 'standard' : uses the standard definition to compute the dynamic modes, `F = U * W`.
Expand All @@ -66,10 +66,10 @@ def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,
Matrix containing the dynamic modes of shape `(m, n-1)` or `(m, k)`.

b : array_like, if `return_amplitudes=True`
1-D array containing the amplitudes of length `min(n-1, k)`.
1-D array containing the amplitudes of length `min(n-1, rank)`.

V : array_like, if `return_vandermonde=True`
Vandermonde matrix of shape `(n-1, n-1)` or `(k, n-1)`.
Vandermonde matrix of shape `(n-1, n-1)` or `(rank, n-1)`.

omega : array_like
Time scaled eigenvalues: `ln(l)/dt`.
Expand Down Expand Up @@ -97,8 +97,8 @@ def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,
raise ValueError('A.dtype must be one of %s, not %s'
% (' '.join(_VALID_DTYPES), A.dtype))

if k is not None and (k < 1 or k > n - 1):
raise ValueError('k must be > 1 and less than n - 1')
if rank is not None and (rank < 1 or rank > n ):
raise ValueError('rank must be > 1 and less than n')

#Split data into lef and right snapshot sequence
X = A[:, :(n-1)] #pointer
Expand All @@ -108,10 +108,10 @@ def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,
U, s, Vh = linalg.svd(X, compute_uv=True, full_matrices=False,
overwrite_a=False, check_finite=True)

if k is not None:
U = U[:, :k]
s = s[:k]
Vh = Vh[:k, :]
if rank is not None:
U = U[:, :rank]
s = s[:rank]
Vh = Vh[:rank, :]

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Solve the LS problem to find estimate for M using the pseudo-inverse
Expand Down Expand Up @@ -156,9 +156,8 @@ def dmd(A, dt=1, k=None, modes='exact', return_amplitudes=False,
return result


def rdmd(A, dt=1, k=None, p=10, l=None, q=2, sdist='uniform', single_pass=False,
return_amplitudes=False, return_vandermonde=False, order=True,
random_state=None):
def rdmd(A, rank, dt=1, oversample=10, n_subspace=2, return_amplitudes=False,
return_vandermonde=False, order=True, random_state=None):
"""Randomized Dynamic Mode Decomposition.

Dynamic Mode Decomposition (DMD) is a data processing algorithm which
Expand All @@ -167,34 +166,29 @@ def rdmd(A, dt=1, k=None, p=10, l=None, q=2, sdist='uniform', single_pass=False,
The modes are ordered corresponding to the amplitudes stored in the diagonal
matrix `B`. `V` is a Vandermonde matrix describing the temporal evolution.

The quality of the approximation can be controlled via the oversampling
parameter `oversample` and `n_subspace` which specifies the number of
subspace iterations.


Parameters
----------
A : array_like
Real/complex input matrix `a` with dimensions `(m, n)`.
A : array_like, shape `(m, n)`.
Input array.

rank : integer
Target rank. Best if `rank << min{m,n}`

dt : scalar or array_like
Factor specifying the time difference between the observations.

k : int
If `k < (n-1)` low-k Dynamic Mode Decomposition is computed.
oversample : integer, optional (default: 10)
Controls the oversampling of column space. Increasing this parameter
may improve numerical accuracy.

p : integer, default: `p=10`.
Parameter to control oversampling of column space.

l : integer, default: `l=2*p`.
Parameter to control oversampling of row space. Only relevant if
single_pass == True.

q : int, optional
Number of subspace iterations to perform. Only relevant if
single_pass == False

sdist : str `{'uniform', 'normal'}`
Specify the distribution of the sensing matrix `S`.

single_pass : bool
If single_pass == True, perfom single pass of algorithm.
n_subspace : integer, default: 2.
Parameter to control number of subspace iterations. Increasing this
parameter may improve numerical accuracy.

return_amplitudes : bool `{True, False}`
True: return amplitudes in addition to dynamic modes.
Expand All @@ -214,13 +208,13 @@ def rdmd(A, dt=1, k=None, p=10, l=None, q=2, sdist='uniform', single_pass=False,
Returns
-------
F : array_like
Matrix containing the dynamic modes of shape `(m, n-1)` or `(m, k)`.
Matrix containing the dynamic modes of shape `(m, rank)`.

b : array_like, if `return_amplitudes=True`
1-D array containing the amplitudes of length `min(n-1, k)`.
1-D array containing the amplitudes of length `min(n-1, rank)`.

V : array_like, if `return_vandermonde=True`
Vandermonde matrix of shape `(n-1, n-1)` or `(k, n-1)`.
Vandermonde matrix of shape `(rank, n-1)`.

omega : array_like
Time scaled eigenvalues: `ln(l)/dt`.
Expand All @@ -233,13 +227,12 @@ def rdmd(A, dt=1, k=None, p=10, l=None, q=2, sdist='uniform', single_pass=False,
(available at `arXiv <https://arxiv.org/abs/1609.00048>`_).
"""
# Compute QB decomposition
Q, B = rqb(A, k=k, p=p, l=l, q=q, sdist=sdist, single_pass=single_pass,
random_state=random_state)
Q, B = rqb(A, rank, oversample=oversample, n_subspace=n_subspace, random_state=random_state)

# only difference is we need to premultiply F from dmd
# vandermonde is basically already computed
# TODO: factor out the rest so no code is repeated
F, V, omega = dmd(B, dt=dt, k=k, modes='standard',return_amplitudes=False,
F, V, omega = dmd(B, rank=rank, dt=dt, modes='standard', return_amplitudes=False,
return_vandermonde=True, order=order)

#Compute DMD Modes
Expand Down
Loading