首先，永远不要计算显式逆，尽管有时你无法避免它。然而，大多数时候，解决方案可以重述为向量与逆矩阵的乘积。在这种情况下，存在几种方法，它们不会明确计算逆。

也就是说，常见的正则化技术是截断奇异值分解（参见截断 svd）或tikhonov 正则化。这两种技术本质上都是过滤小的奇异值。这是必要的，因为小奇异值引入到解向量中的误差可能会在一定程度上占主导地位，即解向量的所有数字都关闭。

现在手头有两种不同的技术，下一个问题是，如何选择正则化参数？这个想法是平衡解向量中的误差与残差的误差。对于我的问题，l 曲线标准通常工作得很好。我强烈推荐那篇论文。它描述了从数据矩阵的 svd 开始计算 l 曲线标准的所有步骤。

我还对截断的 svd 和 tikhonov 正则化进行了快速比较。如果根据 l 曲线标准选择截止/正则化参数，两者都会给出非常相似的结果。稍后我可能会用一个简短的 python 程序更新我的帖子，展示这两种技术。

正如所承诺的那样，一个代码片段计算了一个病态问题的 tikhonov 参数。该tikhonov_lcurve例程使用与论文中相同的符号。

# -*- coding: utf8 -*-
'''
Created on Mar 10, 2013

'''
import numpy as np
import matplotlib.pyplot as pl
import os
import sys

def fun(x):
    # witch of agnesi funcion
    # because of the two complex poles near the x-axes at z = +/- I epsilon
    # many chebyshev modes are necessary for approximation
    eps = 0.03
    return eps ** 2 / (x ** 2 + eps ** 2) 

def main():
    # setting up an ill-posed problem
    # approximate the witch of agnesi function
    # up to order 100 with 100 random points
    np.random.seed(1)
    N = 100           # number of chebyshev modes
    NN = 100          # number of data points which are not gauss quadrature points

    x = 2 * np.random.rand(NN) - 1 # 100 random points in the interval -1,1

    n = np.arange(N)
    nn, xx = np.meshgrid(n, x)
    L = np.cos(nn * np.arccos(xx))
    b = fun(x)

    # ill posed problem stated: L c = b
    # L in this case is my data matrix

    U, sigma, Vh = np.linalg.svd(L, full_matrices=False) # do the svd  
    print 'condition of L:', sigma[0] / sigma[-1]

    beta = np.dot(U.T, b) # setup right side of the ill-posed problem
    x1, y1 = truncation_lcurve(U.copy(), sigma.copy(), Vh.copy(), beta.copy())
    x2, y2, xb, yb = tikhonov_lcurve(U.copy(), sigma.copy(), Vh.copy(), beta.copy())

    # plot the balance between norm of solution vector and norm of residuum vector
    fig = pl.figure()
    ax = fig.add_subplot(111)
    ax.loglog(x2, y2, '.-',label='Tikhonov')
    ax.loglog(x1, y1, '.-',label='Truncated SVD')
    ax.loglog(xb, yb, 'o',label='Best Tikhonov Parameter')
    ax.set_xlim(1e-8, 1)
    ax.set_ylim(0.01, 1e10)
    ax.set_xlabel('|Ax-b|')
    ax.set_ylabel('|x|')
    pl.legend()
    pl.show()

def truncation_lcurve(U, sigma, Vh, beta):
    l = np.arange(90) + 10  # truncate singular values from 100 down to 10 and check the effect on the lcurve

    ll, ss = np.meshgrid(l, sigma)
    f = (0 * ss + 1)
    for i, v in enumerate(l):
        s = np.ones(sigma.shape)
        s[v:] = 0.                # set singular values after cut off to zero == truncation == singular value filter
        f[:, i] = s

    _, bbeta = np.meshgrid(l, beta)

    # xl contains the regularized solution for all l trunctions
    xl = np.dot(Vh.T, f * bbeta / ss)


    # the following calculations are necessary for the l-curve criterion, see publication
    nxl = np.sum(np.abs(xl) ** 2, axis=0) ** (1. / 2)# norm of the solution vector
    axmb = np.sum(np.abs((1 - f) * bbeta) ** 2, axis=0) ** (1. / 2)# norm of the resdiuum vector

    return axmb, nxl


def tikhonov_lcurve(U, sigma, Vh, beta):
    # as shown in the publication, the kink is between the lowest and highest singular value
    # let's try 1000 regularization parameters in the interval of the singular values
    l = np.logspace(np.log10(sigma[-1]) - 1, np.log10(sigma[0]) + 1, 1000) 

    ll, ss = np.meshgrid(l, sigma)
    # regularize the singular value with ll**2, note if ss >> ll it has no effect, 
    # if ss is of order of ll however, damping kicks in
    # that's why one calls both regularization technique
    # a filter for singular values
    f = ss ** 2 / (ss ** 2 + ll ** 2)
    _, bbeta = np.meshgrid(l, beta)

    # xl contains the regularized solution for all l parameters
    xl = np.dot(Vh.T, f * bbeta / ss)

    # the following calculations are necessary for the l-curve criterion, see publication
    # e.g. kappa is the curvate of the l-curve on the double log scale
    nxl = np.sum(np.abs(xl) ** 2, axis=0) ** (1. / 2)  # norm of the solution vector
    axmb = np.sum(np.abs((1 - f) * bbeta) ** 2, axis=0) ** (1. / 2) # norm of the resdiuum vector

    eta = nxl ** 2
    rho = axmb ** 2

    eta_p = -4 / l * np.sum((1 - f) * f ** 2 * bbeta ** 2 / ss ** 2, axis=0)
    kappa = 2 * eta * rho / eta_p * (l ** 2 * eta_p * rho + 2 * l * eta * rho + l ** 4 * eta * eta_p) / (l ** 2 * eta ** 2 + rho ** 2) ** (3. / 2)
    bi = kappa.argmin()
    lb = l[bi]
    print 'Best Tikhonov Parameter ', lb

    # xbl contains the solution where the tikhonov regularization
    # is near the kink of the l-curve
    xbl = np.dot(Vh.T, sigma / (sigma ** 2 + lb ** 2) * beta)
    nxbl = np.linalg.norm(xbl)
    axmb_xbl = np.linalg.norm((1 - sigma ** 2 / (sigma ** 2 + lb ** 2)) * beta)

    return axmb, nxl, axmb_xbl, nxbl

if __name__ == '__main__':
    main()