Merge pull request #96 from JuliaImageRecon/nh/tikhMat

Tikhonov matrix bug fixes for Kaczmarz extended system of equations

src/Kaczmarz.jl

@@ -155,7 +155,7 @@ function init!(solver::Kaczmarz, state::KaczmarzState{T, vecT}, b::vecT; x0 = 0)
 
   state.u .= b
   if λ_ isa AbstractVector
-    state.ɛw = 0
+    state.ɛw = one(T)
   else
     state.ɛw = sqrt(λ_)
   end
@@ -240,12 +240,18 @@ function initkaczmarz(A,λ)
   return A, denom, rowindex
 end
 function initkaczmarz(A, λ::AbstractVector)
+  # Instead of ||Ax - b||² + λ||x||² we solve ||Ax - b||² + λ||Lx||² where L is a diagonal matrix
+  # See Per Christian Hansen: Rank-Deficient and Discrete Ill-Posed Problems, Chapter 2.3 Transformation to Standard Form
+  # -> ||Âc - u||² + λ||c||² with  = AL⁻¹, c = Lx
+  # We put this into the standard extended system of equation with λ = 1
+  # In the end we need to multiply the solution with L⁻¹ 
   λ = real(eltype(A)).(λ)
   A = initikhonov(A, λ)
-  return initkaczmarz(A, 0)
+  return initkaczmarz(A, one(eltype(λ)))
 end
 
-initikhonov(A, λ) = transpose((1 ./ sqrt.(λ)) .* transpose(A)) # optimize structure for row access
+# A * inv(λ), specialised for diagm(λ)
+initikhonov(A, λ::AbstractVector) =  transpose((1 ./ sqrt.(λ)) .* transpose(A)) # optimize structure for row access
 initikhonov(prod::ProdOp{Tc, <:WeightingOp, matT}, λ) where {T, Tc<:Union{T, Complex{T}}, matT} = ProdOp(prod.A, initikhonov(prod.B, λ))
 ### kaczmarz_update! ###
 

test/testKaczmarz.jl

@@ -37,18 +37,6 @@ Random.seed!(12345)
 
       # Test Tikhonov regularization matrix
       @testset "Kaczmarz Tikhonov matrix" begin
-        A = rand(3, 2) + im * rand(3, 2)
-        x = rand(2) + im * rand(2)
-        b = A * x
-
-        regMatrix = rand(2) # Tikhonov matrix
-
-        solver = Kaczmarz
-        S = createLinearSolver(solver, arrayType(A), iterations=200, reg=[L2Regularization(arrayType(regMatrix))])
-        x_approx = Array(solve!(S, arrayType(b)))
-        #@info "Testing solver $solver ...: $x  == $x_approx"
-        @test norm(x - x_approx) / norm(x) ≈ 0 atol = 0.1
-
         ## Test spatial regularization
         M = 12
         N = 8
@@ -62,17 +50,25 @@ Random.seed!(12345)
 
         # @show A, x, regMatrix
         # use regularization matrix
-
-        S = createLinearSolver(solver, arrayType(A), iterations=100, reg=[L2Regularization(arrayType(regMatrix))])
+        S = createLinearSolver(Kaczmarz, arrayType(A), iterations=100, reg=[L2Regularization(arrayType(regMatrix))])
         x_matrix = Array(solve!(S, arrayType(b)))
 
         # use standard reconstruction
-        S = createLinearSolver(solver, arrayType(A * Diagonal(1 ./ sqrt.(regMatrix))), iterations=100)
+        S = createLinearSolver(Kaczmarz, arrayType(A * Diagonal(1 ./ sqrt.(regMatrix))), reg = [L2Regularization(1.0)], iterations=100)
         x_approx = Array(solve!(S, arrayType(b))) ./ sqrt.(regMatrix)
 
         # test
         #@info "Testing solver $solver ...: $x_matrix  == $x_approx"
         @test norm(x_approx - x_matrix) / norm(x_approx) ≈ 0 atol = 0.1
+
+        # Compare reg. matrix of equal elements to standard reco
+        λ = rand()
+        S = createLinearSolver(Kaczmarz, arrayType(A), iterations=100, reg=[L2Regularization(λ)])
+        x_standard = Array(solve!(S, arrayType(b)))
+
+        S = createLinearSolver(Kaczmarz, arrayType(A), iterations=100, reg=[L2Regularization(arrayType(fill(λ, N)))])
+        x_matrix = Array(solve!(S, arrayType(b)))
+        @test isapprox(x_standard, x_matrix)
       end
 
       @testset "Kaczmarz Weighting Matrix" begin