In [1]:
import numpy as np
from dec.schur import schur, schurλ, bmat, rand
from dec.spectral import I_diag
In [2]:
n, m = 5, 2
B = np.arange(n*m).reshape((n,m))
C = np.arange(n*m).reshape((m,n))
D = np.arange(m*m).reshape((m,m))
print(B)
print(C)
print(D)
print(B.dot(np.eye(B.shape[-1], B.shape[-1])))
print(C.dot(np.eye(C.shape[-1], C.shape[-1])))
In [3]:
def fourier_sin(N, a, b):
r'''
Corresponds to :math:`f(x) \mapsto \int_{x+a}^{x+b} f(\xi) \sin(\xi) d\xi`
'''
I = I_diag(N+2, a, b) / 2j
def K(x):
assert x.shape == (N,)
x = np.array(x, dtype=np.complex)
x = np.hstack([[0], x, [0]])
x = (np.roll(x,+1) - np.roll(x,-1))
x *= I
rslt = x[1:-1]
rslt[ 0] += x[-1]
rslt[-1] += x[0]
return rslt
def Kinv(x):
# Make sure type is coerced to complex, otherwise numpy ignores the complex parts
# and reverts to reals.
assert x.shape == (N,)
x = np.array(x.copy(), dtype=np.complex)
x /= I[1:-1]
if (np.isclose(I[0], I[N]) or
np.isclose(I[1], I[N+1]) or
np.isclose(I[0]*I[1], I[N]*I[N+1])):
raise ValueError("Singular operator.")
y = np.zeros(N, dtype=np.complex)
# The computations below are essentially Schur's complement?
E = np.sum(x[::2]); O = np.sum(x[1::2])
if N % 2 == 0:
y[0] = O/(1-I[0]/I[N])
y[-1] = E/(I[N+1]/I[1]-1)
else:
y[0] = (I[1]/I[N+1]*E+O)/(1-I[1]*I[0]/I[N]/I[N+1])
y[-1] = (I[N]/I[0]*E+O)/(I[N]*I[N+1]/I[0]/I[1]-1)
x[0] -= y[-1]*I[N+1]/I[1]
x[-1] -= -y[0]*I[0]/I[N]
# This should be the crux of the inverse
x = np.hstack([[-y[0]], x , [y[-1]]])
y[::2] = -np.cumsum(x[::2])[:-1]
y[1::2] = np.cumsum(x[1::2][::-1])[:-1][::-1]
return y
return K, Kinv
In [4]:
N = 4
K, Kinv = fourier_sin(N, 0, 0.1)
x = np.ones(N)
y = K(x)
assert np.allclose(x, Kinv(y))
In [5]:
import dec.spectral as sp
from dec.helper import to_matrix
def A_bidiagonal(x):
f = np.hstack([ [-x[1]], x[:-2]-x[2:], [x[-2]] ])
return f
to_matrix(A_bidiagonal, 6)
#print(np.linalg.inv(to_matrix(Aa, 6)))
Out[5]:
In [6]:
def fourier_sin_schur_old(N, a, b):
I = sp.I_diag(N+2, a, b) / 2j
zeros = lambda n: np.zeros(n, dtype=np.complex)
n = N - 2
m = 2
def A(x):
return A_bidiagonal(x)*I[2:-2]
def Ainv(f):
raise NotImplemented
def B(λ):
f = zeros(n)
f[0] = I[ 2]*λ[0]
f[-1] = -I[-3]*λ[1]
return f
def C(x):
g = zeros(m)
g[0] = -I[ 1]*x[ 0]
g[1] = I[-2]*x[-1]
return g
def D(λ):
g = zeros(m)
g[0] = I[-1]*λ[1]
g[1] = -I[ 0]*λ[0]
return g
def K(x):
assert len(x) == N
x = np.array(x, dtype=np.complex)
x, λ = x[1:-1], x[[0, -1]]
f = A(x) + B(λ)
g = C(x) + D(λ)
return np.hstack([ [g[0]], f, [g[1]] ])
Minv = None
def Kinv(f):
raise NotImplemented
return K, Kinv
for N in range(4, 9):
K_, Kinv_ = fourier_sin(N, 0, 1)
K, Kinv = fourier_sin_schur_old(N, 0, 1)
x = np.random.rand(N)
assert np.allclose(K(x), K_(x))
In [7]:
def A_lower_triangular(x):
if len(x)<3: return x
return np.hstack([ [x[0], x[1]], x[2:]-x[:-2] ])
to_matrix(A_lower_triangular, 5)
Out[7]:
In [8]:
np.linalg.inv( to_matrix(A_lower_triangular, 5) )
Out[8]:
In [9]:
def A_lower_triangular_inv(x):
if len(x)<3: return x
y = np.zeros_like(x)
y[::2] = np.cumsum(x[::2])
y[1::2] = np.cumsum(x[1::2])
return y
to_matrix(A_lower_triangular_inv, 5)
Out[9]:
In [10]:
for n in range(4, 10):
x = np.random.rand(n)
assert np.allclose(x, A_lower_triangular(A_lower_triangular_inv(x)))
assert np.allclose(x, A_lower_triangular_inv(A_lower_triangular(x)))
In [11]:
def fourier_sin_schur(N, a, b):
I = sp.I_diag(N+2, a, b) / 2j
zeros = lambda n: np.zeros(n, dtype=np.complex)
n = N - 2
m = 2
def A(x):
return A_lower_triangular(-x)*I[2:-2]
def Ainv(f):
return A_lower_triangular_inv(-f/I[2:-2])
def B(λ):
f = zeros(n)
f[0] = I[2]*λ[0]
f[1] = I[3]*λ[1]
return f
def C(x):
g = zeros(m)
g[0] = I[-1]*x[-1]
g[1] = I[-2]*x[-2]
return g
def D(λ):
g = zeros(m)
g[0] = -I[1]*λ[1]
g[1] = -I[0]*λ[0]
return g
def K(x):
assert len(x) == N
x = np.array(x, dtype=np.complex)
x, λ = x[2:], x[:2]
f = A(x) + B(λ)
g = C(x) + D(λ)
return np.hstack([ [g[0]], f, [g[1]] ])
Minv = schurλ(Ainv, B, C, D, n, m)
def Kinv(f):
assert len(f) == N
f, g = f[1:-1], f[[0,-1]]
x, λ = Minv(f, g)
return np.hstack([λ, x])
return K, Kinv
for N in range(4, 9):
K_, Kinv_ = fourier_sin(N, 0, 1)
K, Kinv = fourier_sin_schur(N, 0, 1)
x = np.random.rand(N)
assert np.allclose(K(x), K_(x))
assert np.allclose(Kinv(x), Kinv_(x))
In [12]:
def fourier_sin_schur_simple(N, a, b):
I = sp.I_diag(N+2, a, b) / 2j
zeros = lambda n: np.zeros(n, dtype=np.complex)
n = N - 2
m = 2
K, _ = fourier_sin(N, a, b)
zn = lambda: zeros(n)
zm = lambda: zeros(m)
def Kb(x=zn(), λ=zm()):
y = K(np.hstack([λ, x]))
f, g = y[1:-1], y[[0, -1]]
return f, g
def Ainv(f):
return A_lower_triangular_inv(-f/I[2:-2])
def A(x): return Kb(x=x)[0]
def B(λ): return Kb(λ=λ)[0]
def C(x): return Kb(x=x)[1]
def D(λ): return Kb(λ=λ)[1]
Minv = schurλ(Ainv, B, C, D, n, m)
def Kinv(f):
assert len(f) == N
f, g = f[1:-1], f[[0,-1]]
x, λ = Minv(f, g)
return np.hstack([λ, x])
return K, Kinv
for N in range(4, 9):
K_, Kinv_ = fourier_sin(N, 0, 1)
K, Kinv = fourier_sin_schur_simple(N, 0, 1)
x = np.random.rand(N)
assert np.allclose(K(x), K_(x))
assert np.allclose(Kinv(x), Kinv_(x))
In [13]:
def I_diag_sy(N, a, b):
from sympy import symbols
l = symbols(['I{}'.format(i) for i in range(N)])
return np.array(l)
#def I_diag_sy(N, a, b):
# from sympy import I, exp, Integer, symbols
# l = []
# for n in map(int, sp.freq(N)):
# if n == 0:
# l.append(b-a)
# else:
# n = Integer(n)
# l.append( (exp(I*n*b) - exp(I*n*a))/(I*n) )/2/I
# return np.array(l)
def fourier_sin_schur_sym(a, b):
from sympy import I, exp, Integer
def K(x):
N = x.shape[0]
x = np.hstack([[0], x, [0]])
x = (np.roll(x,+1) - np.roll(x,-1))
x *= I_diag_sy(N+2, a, b)
rslt = x[1:-1]
rslt[ 0] += x[-1]
rslt[-1] += x[0]
return rslt
return K, None
import sympy as sy
a, b = sy.symbols('a b')
K, Kinv = fourier_sin_schur_sym(a, b)
In [14]:
X = np.array(sy.symbols('λ0 x0 x1 x2 x3 x4 x5 x6 x7 x8 λ1'))
K(X)
Out[14]:
In [15]:
# This is better because we have a lower triangular matrix.
X = np.array(sy.symbols('λ0 λ1 x0 x1 x2 x3 x4 x5 x6 x7 x8'))
K(X)
Out[15]:
In [16]:
X = np.array(sy.symbols('λ0 λ1'))
K(X)
Out[16]:
In [17]:
X = np.array(sy.symbols('λ0 λ1 x0'))
K(X)
Out[17]:
In [18]:
A_lower_triangular(I_diag_sy(5, a, b))
Out[18]:
In [19]:
N = 1000
K, Kinv = fourier_sin(N, 0, 0.1)
K_, Kinv_ = fourier_sin_schur(N, 0, 0.1)
K__, Kinv__ = fourier_sin_schur_simple(N, 0, 0.1)
np.random.seed(1)
x = np.random.rand(N)
assert np.allclose(K_(x), K(x))
assert np.allclose(K__(x), K(x))
assert np.allclose(Kinv_(x), Kinv(x))
assert np.allclose(Kinv__(x), Kinv(x))
%timeit K(x)
%timeit Kinv(x)
%timeit K_(x)
%timeit Kinv_(x)
%timeit K__(x)
%timeit Kinv__(x)
In [20]: