5
$\begingroup$

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:

  1. Pick 2000
  2. calculate the weights via qr
  3. 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()                                                                 
$\endgroup$
  • $\begingroup$ On your steps (1) and (2), do you mean you are using rank revealing QR? $\endgroup$ – Reid.Atcheson Apr 2 '13 at 7:38
  • $\begingroup$ @Reid.Atcheson No I don't. In my toy system I do a svd. I can access the rank via np.sum(s>s[0]*eps) where eps is the numerical cut off parameter. I see, that the numerical rank deficiency is different from the number of negative weights. Does this mean I replace too many? $\endgroup$ – Bort Apr 2 '13 at 13:15
  • $\begingroup$ I was going to say, a trick I use a lot when brute force isn't too slow is just to sample something to death and then do a rank revealing QR (nice trick in finding good quadrature points on weird shapes). Having the SVD available though is probably better. $\endgroup$ – Reid.Atcheson Apr 3 '13 at 2:38
  • $\begingroup$ Have you already seen the papers by Golub, Elhay, and Kautsky for finding quadrature weights? $\endgroup$ – J. M. Apr 18 '13 at 12:00
  • $\begingroup$ Thanks for pointing me in the right direction. Currently I convinced the fe guy to add grid points add quadrature points which incredibly stabilizes everything. I still look into this and will post a solution on my toy system. This might come handy one day. $\endgroup$ – Bort Apr 19 '13 at 21:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.