目前，我在光谱的基础上拟合 PDE 的有限元解决方案。矩阵（$R^{25000\times 2000}$) 的相应线性方程组是高度病态的 ($\kappa \approx 10^{15}$）。但是，我能够找到使用 Tikhonov 正则化的解决方案，该解决方案在域中的其他点上也有少量残留。

现在，在玩玩具系统时，我想出了在所有点的给定子集上找到正交权重的想法，这样我就必须从 100000 FE 中找到 2000 点，而不是过度拟合 25000 点点，使得它们具有正正交权重。

您建议从 100000 中选择这 2000 点的策略是什么？

目前我执行以下操作：

1. 选择 2000
2. 通过 qr 计算权重
3. 选择所有具有负权重的点并用新点替换它们

在我的具有切比雪夫多项式的一维玩具系统上，我需要为 25 种模式的频谱基础绘制大约 100000 个新点，直到我的所有权重都是正数。所以这种方法不会奏效。

非常欢迎任何建议或文献链接。谢谢。

这是我的玩具系统：

# 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()