From 63b2f778d28729fcf533a2fedadee258a599884b Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 26 Nov 2024 13:29:46 +0100 Subject: [PATCH 1/4] Fix transformation in Kaczmarz with tikhonov matrix --- src/Kaczmarz.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Kaczmarz.jl b/src/Kaczmarz.jl index 3584c85..882c1b0 100644 --- a/src/Kaczmarz.jl +++ b/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! ### From 6e5c1dbf446c01df4e9e31ddf111cad1496714fd Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 26 Nov 2024 13:30:04 +0100 Subject: [PATCH 2/4] Update Kaczmarz tikhonov tests --- test/testKaczmarz.jl | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/test/testKaczmarz.jl b/test/testKaczmarz.jl index 6839583..8e9f256 100644 --- a/test/testKaczmarz.jl +++ b/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))]) x_matrix = Array(solve!(S, arrayType(b))) # use standard reconstruction - S = createLinearSolver(solver, arrayType(A * Diagonal(1 ./ sqrt.(regMatrix))), iterations=100) + S = createLinearSolver(solver, 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(solver, arrayType(A), iterations=100, reg=[L2Regularization(λ)]) + x_standard = Array(solve!(S, arrayType(b))) + + S = createLinearSolver(solver, arrayType(A), iterations=100, reg=[L2Regularization(fill(λ, N))]) + x_matrix = Array(solve!(S, arrayType(b))) + @test isapprox(x_standard, x_matrix) end @testset "Kaczmarz Weighting Matrix" begin From 178a58190862cb0f13ad42fb0debdd25054d6e2a Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 26 Nov 2024 15:18:20 +0100 Subject: [PATCH 3/4] Fix Kaczmarz tikhonov test --- test/testKaczmarz.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/testKaczmarz.jl b/test/testKaczmarz.jl index 8e9f256..a877069 100644 --- a/test/testKaczmarz.jl +++ b/test/testKaczmarz.jl @@ -50,11 +50,11 @@ 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))), reg = [L2Regularization(1.0)], 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 @@ -63,10 +63,10 @@ Random.seed!(12345) # Compare reg. matrix of equal elements to standard reco λ = rand() - S = createLinearSolver(solver, arrayType(A), iterations=100, reg=[L2Regularization(λ)]) + S = createLinearSolver(Kaczmarz, arrayType(A), iterations=100, reg=[L2Regularization(λ)]) x_standard = Array(solve!(S, arrayType(b))) - S = createLinearSolver(solver, arrayType(A), iterations=100, reg=[L2Regularization(fill(λ, N))]) + S = createLinearSolver(Kaczmarz, arrayType(A), iterations=100, reg=[L2Regularization(fill(λ, N))]) x_matrix = Array(solve!(S, arrayType(b))) @test isapprox(x_standard, x_matrix) end From d5ccb8612b4c7167aa4f54d85b42b53ee73a7b97 Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 26 Nov 2024 16:11:30 +0100 Subject: [PATCH 4/4] Fix Kaczmarz tikhonov test for GPUs --- test/testKaczmarz.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testKaczmarz.jl b/test/testKaczmarz.jl index a877069..e4dd3ee 100644 --- a/test/testKaczmarz.jl +++ b/test/testKaczmarz.jl @@ -66,7 +66,7 @@ Random.seed!(12345) 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(fill(λ, N))]) + 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