# SVD regularization - ray 2D tomography

Sunny day today, isn't it? Please, I need help with my problem. I have written a program to do 2D ray tomography, according to this paper.

For the result, I use formula (4.15) from the paper. Now I need to do SVD regularization on the result, in order to get good results for lower count of rays (lets say I will have $N$ boxes and $N$ rays). Now it is purely impossible with the formula (4.15)

I am not sure how to understand the regularization - I think it should be only to add something to the diagonal, right? Please can you help me with it?

I have found the paper "Regularized matrix computation" from AE Yagle, but I'm not sure how to implement it in my case. Many thanks.

• It's not sunny here! :-) – Michael Grant May 21 '13 at 15:21

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.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()
• "First, never calculate an explicit inverse, allthough sometimes you cannot avoid it." It is true that formulas that are written in terms of matrix inverses can often be cleverly computed without actually inverting a matrix, and the folklore in scientific computing is that it is desirable to be clever in this way. On the other hand, someone has posted the following heretical paper onto arxiv claiming that avoiding inverses is less important than our culture would suggest arxiv.org/abs/1201.6035 – clipper May 24 '13 at 17:32