Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63727257
pod_engine.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, May 22, 02:52
Size
4 KB
Mime Type
text/x-python
Expires
Fri, May 24, 02:52 (2 d)
Engine
blob
Format
Raw Data
Handle
17811546
Attached To
R6746 RationalROMPy
pod_engine.py
View Options
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import
numpy
as
np
from
copy
import
copy
from
rrompy.utilities.types
import
Np1D
,
Np2D
,
Tuple
__all__
=
[
'PODEngine'
]
class
PODEngine
:
"""
POD engine for general matrix orthogonalization.
Args/Attributes:
M: Square matrix associated to scalar product.
"""
def
__init__
(
self
,
M
:
Np2D
):
self
.
M
=
M
def
norm
(
self
,
a
:
Np1D
)
->
float
:
"""Compute norm of a vector."""
return
np
.
sqrt
(
np
.
abs
(
a
.
conj
()
.
T
.
dot
(
self
.
M
.
dot
(
a
))))
def
GS
(
self
,
a
:
Np1D
,
Q
:
Np2D
,
n
:
int
=
None
)
->
Tuple
[
Np1D
,
Np1D
]:
"""
Compute 1 Gram-Schmidt step with given projector.
Args:
a: vector to be projected;
Q: orthogonal projection matrix;
n: number of columns of Q to be considered.
Returns:
Resulting normalized vector, coefficients of a wrt the updated
basis.
"""
if
n
is
None
:
n
=
Q
.
shape
[
1
]
Q
=
Q
[:,
:
n
]
r
=
np
.
zeros
((
n
+
1
,),
dtype
=
a
.
dtype
)
if
n
>
0
:
for
j
in
range
(
2
):
# twice is enough!
nu
=
Q
.
conj
()
.
T
.
dot
(
self
.
M
.
dot
(
a
))
a
=
a
-
Q
.
dot
(
nu
)
r
[:
-
1
]
=
r
[:
-
1
]
+
nu
r
[
-
1
]
=
self
.
norm
(
a
)
if
np
.
isclose
(
np
.
abs
(
r
[
-
1
]),
0.
):
r
[
-
1
]
=
1.
a
=
a
/
r
[
-
1
]
return
a
,
r
def
QRGramSchmidt
(
self
,
A
:
Np2D
,
only_R
:
bool
=
False
)
->
Tuple
[
Np1D
,
Np1D
]:
"""
Compute QR decomposition of a matrix through Gram-Schmidt method.
Args:
A: matrix to be decomposed;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting orthogonal and upper-triangular factors.
"""
N
=
A
.
shape
[
1
]
Q
=
np
.
zeros_like
(
A
,
dtype
=
A
.
dtype
)
R
=
np
.
zeros
((
N
,
N
),
dtype
=
A
.
dtype
)
for
k
in
range
(
N
):
Q
[:,
k
],
R
[:
k
+
1
,
k
]
=
self
.
GS
(
A
[:,
k
],
Q
,
k
)
if
only_R
:
return
R
return
Q
,
R
def
QRHouseholder
(
self
,
A
:
Np2D
,
Q0
:
Np2D
=
None
,
only_R
:
bool
=
False
)
->
Tuple
[
Np1D
,
Np1D
]:
"""
Compute QR decomposition of a matrix through Householder method.
Args:
A: matrix to be decomposed;
Q0(optional): initial orthogonal guess for Q; defaults to random;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting (orthogonal and )upper-triangular factor(s).
"""
B
=
copy
(
A
)
N
=
B
.
shape
[
1
]
V
=
np
.
zeros_like
(
B
,
dtype
=
B
.
dtype
)
R
=
np
.
zeros
((
N
,
N
),
dtype
=
B
.
dtype
)
if
Q0
is
None
:
Q
=
np
.
zeros_like
(
B
,
dtype
=
B
.
dtype
)
+
np
.
random
.
randn
(
*
(
B
.
shape
))
else
:
Q
=
copy
(
Q0
)
for
k
in
range
(
N
):
if
Q0
is
None
:
Q
[:,
k
],
_
=
self
.
GS
(
Q
[:,
k
],
Q
,
k
)
a
=
B
[:,
k
]
R
[
k
,
k
]
=
self
.
norm
(
a
)
alpha
=
Q
[:,
k
]
.
conj
()
.
T
.
dot
(
self
.
M
.
dot
(
a
))
if
np
.
isclose
(
np
.
abs
(
alpha
),
0.
):
s
=
1.
else
:
s
=
-
alpha
/
np
.
abs
(
alpha
)
Q
[:,
k
]
=
s
*
Q
[:,
k
]
V
[:,
k
],
_
=
self
.
GS
(
R
[
k
,
k
]
*
Q
[:,
k
]
-
a
,
Q
,
k
)
J
=
np
.
arange
(
k
+
1
,
N
)
vtB
=
V
[:,
k
]
.
conj
()
.
dot
(
self
.
M
.
dot
(
B
[:,
J
]))
B
[:,
J
]
=
B
[:,
J
]
-
2
*
np
.
outer
(
V
[:,
k
],
vtB
)
R
[
k
,
J
]
=
Q
[:,
k
]
.
conj
()
.
T
.
dot
(
self
.
M
.
dot
(
B
[:,
J
]))
B
[:,
J
]
=
B
[:,
J
]
-
np
.
outer
(
Q
[:,
k
],
R
[
k
,
J
])
if
only_R
:
return
R
for
k
in
range
(
N
-
1
,
-
1
,
-
1
):
J
=
np
.
arange
(
k
,
N
)
vtQ
=
V
[:,
k
]
.
conj
()
.
T
.
dot
(
self
.
M
.
dot
(
Q
[:,
J
]))
Q
[:,
J
]
=
Q
[:,
J
]
-
2
*
np
.
outer
(
V
[:,
k
],
vtQ
)
return
Q
,
R
Event Timeline
Log In to Comment