Currently, I fit a Finite Element solution of a PDE on a spectral basis. The matrices ($R^{25000\times 2000}$) of the corresponding system of linear equations are highly ill-conditioned ($\kappa \approx 10^{15}$). However, I am able to find a solution with Tikhonov-regularization, which has also a small residuum on other points in the domain.
Now, playing around with a toy system, I came up with the idea to find the quadrature weights on a given subset of all points, such that instead of over-fitting with 25000 points, I have to find 2000 points out of my 100000 FE points such that they have positive quadrature weights.
What strategy do you suggest for selecting those 2000 points from the 100000?
Currently I do the following:
- Pick 2000
- calculate the weights via qr
- select all points with negative weights and replace them with new points
On my 1D toy system with chebyshev polynoms, I need to draw about 100000 new points for a spectral basis of 25 modes till all my weights are positive. So this approach is not going to work well.
Any suggestions or literature links are highly welcome. Thanks.
Here is my toy system:
# encoding: utf8
import numpy as np
import matplotlib.pyplot as pl
def fun(x):
x = np.asarray(x)
y = np.zeros(x.shape)
y[np.abs(x) > 0] = np.exp((-(1. / x[np.abs(x) > 0]) ** 2))
return y
def gauss_extrema(n):
i = np.arange(n)
return np.cos(((2 * (i + 1)) - 1) / (2.*n) * np.pi)
def main():
Nx = 25
nx = Nx
x = 2 * np.random.rand(nx) - 1
# x = gauss_extrema(nx)
m = np.arange(Nx)
mm, xx = np.meshgrid(m, x)
L = np.cos(mm * np.arccos(xx))
sp = np.zeros(nx)
sp[0] = np.pi
w = np.empty(nx)
w.fill(-1)
eps = 0.
i = 1
while np.any(w < eps):
u, s, vh = np.linalg.svd(L.T, full_matrices=False)
if i % 10000 == 0:
print i, s[0] / s[-1], np.sum(w < eps)
beta = np.dot(u.T, sp)
w = np.dot(vh.T, 1. / s * beta)
x[w < eps] = 2 * np.random.rand(np.sum(w < eps)) - 1
mm, xx = np.meshgrid(m, x)
L = np.cos(mm * np.arccos(xx))
i += 1
print w
print x
b = fun(x) # setup right side b
L = np.multiply(L, w) # setup weights on all points on L
u = np.linalg.solve(L, b) # solve L u = b for coefficient vector u
nx = 100
x = np.linspace(-1, 1, nx)
mm, xx = np.meshgrid(m, x)
L = np.cos(mm * np.arccos(xx))
L = np.multiply(L, w)
b = fun(x)
fig = pl.figure()
ax = fig.add_subplot(121)
ax.plot(x, b)
ax.plot(x, np.dot(L, u))
ax = fig.add_subplot(122)
ax.plot(x, np.dot(L, u) - b)
pl.show()
if __name__ == '__main__':
main()
