First, never calculate an explicit inverse, allthough sometimes you cannot avoid it. Most of the time however, the solution can be restated as multiplication of a vector with an inverse matrix. In that case several methods exist, which don't exlpicitely calculate the inverse.
That said, common regularization techniques are truncated singular value decomposition(see truncated svd) or tikhonov regularization. Both techniques essentially filter small singular values.
This is necessary because the error introduced by small singular values to the solution vector can become dominant to that extent, that all digits of the solution vector are off.
Now having two different techniques at hand, the next question is, how to choose the regularization parameter? The idea is to balance the error in the solution vector with the error of the residuum.
For my problems, the l-curve criterion worked usually fine.
I can highly recommend that paper. It describes all steps for a calculation of the l-curve criterion starting from the svd of the data matrix.
I also did a quick comparision between the truncated svd and tikhonov regularization. Both give very similiar results if the cut off / regularization parameters are choosen according to the l-curve criterion. I will probably update my post later with a short python program, demonstrating both techniques.
As promised a code snippet calculating tikhonov parameters of an ill posed problem. The tikhonov_lcurve routine uses the same notation as in the paper.
# -*- 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()
