diff --git a/Project.toml b/Project.toml
index e33269f1..c635eea1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ julia = "1.9"
 StatsBase = "0.33, 0.34"
 FLoops = "0.2"
 VectorizationBase = "0.19, 0.21"
-LinearOperatorCollection = "1.0"
+LinearOperatorCollection = "1.2"
 LinearOperators = "2.3.3"
 FFTW = "1.0"
 
diff --git a/src/CGNR.jl b/src/CGNR.jl
index 69c2ebae..7c48c1eb 100644
--- a/src/CGNR.jl
+++ b/src/CGNR.jl
@@ -12,7 +12,6 @@ mutable struct CGNR{matT,opT,vecT,T,R,PR} <: AbstractKrylovSolver
   αl::T
   βl::T
   ζl::T
-  weights::vecT
   iterations::Int64
   relTol::Float64
   z0::Float64
@@ -46,7 +45,6 @@ function CGNR(A
             ; AHA = A'*A
             , reg = L2Regularization(zero(real(eltype(AHA))))
             , normalizeReg::AbstractRegularizationNormalization = NoNormalization()
-            , weights::AbstractVector = similar(AHA, 0)
             , iterations::Int = 10
             , relTol::Real = eps(real(eltype(AHA)))
             )
@@ -83,7 +81,7 @@ function CGNR(A
 
 
   return CGNR(A, AHA,
-    L2, other, x, x₀, pl, vl, αl, βl, ζl, weights, iterations, relTol, 0.0, normalizeReg)
+    L2, other, x, x₀, pl, vl, αl, βl, ζl, iterations, relTol, 0.0, normalizeReg)
 end
 
 """
@@ -106,16 +104,7 @@ function init!(solver::CGNR, b; x0 = 0)
   end
 
   #x₀ = Aᶜ*rl, where ᶜ denotes complex conjugation
-  if solver.A === nothing
-    !isempty(solver.weights) && @info "weights are being ignored if the backprojection is pre-computed"
-    solver.x₀ .= b
-  else
-    if isempty(solver.weights)
-      mul!(solver.x₀, adjoint(solver.A), b)
-    else
-      mul!(solver.x₀, adjoint(solver.A), b .* solver.weights)
-    end
-  end
+  initCGNR(solver.x₀, solver.A, b)
 
   solver.z0 = norm(solver.x₀)
   copyto!(solver.pl, solver.x₀)
@@ -124,6 +113,10 @@ function init!(solver::CGNR, b; x0 = 0)
   solver.L2 = normalize(solver, solver.normalizeReg, solver.L2, solver.A, b)
 end
 
+initCGNR(x₀, A, b) = mul!(x₀, adjoint(A), b)
+initCGNR(x₀, prod::ProdOp{T, <:WeightingOp, matT}, b) where {T, matT} = mul!(x₀, adjoint(prod.B), b .* prod.A.weights)
+initCGNR(x₀, ::Nothing, b) = x₀ .= b
+
 solverconvergence(solver::CGNR) = (; :residual => norm(solver.x₀))
 
 """
diff --git a/src/Kaczmarz.jl b/src/Kaczmarz.jl
index c751e41f..0a46e5fb 100644
--- a/src/Kaczmarz.jl
+++ b/src/Kaczmarz.jl
@@ -1,7 +1,7 @@
 export kaczmarz
 export Kaczmarz
 
-mutable struct Kaczmarz{matT,T,U,R,RN} <: AbstractRowActionSolver
+mutable struct Kaczmarz{matT,R,T,U,RN} <: AbstractRowActionSolver
   A::matT
   u::Vector{T}
   L2::R
@@ -11,17 +11,15 @@ mutable struct Kaczmarz{matT,T,U,R,RN} <: AbstractRowActionSolver
   rowIndexCycle::Vector{Int64}
   x::Vector{T}
   vl::Vector{T}
-  εw::Vector{T}
+  εw::T
   τl::T
   αl::T
-  weights::Vector{U}
   randomized::Bool
   subMatrixSize::Int64
   probabilities::Vector{U}
   shuffleRows::Bool
   seed::Int64
   iterations::Int64
-  regMatrix::Union{Nothing,Vector{U}} # Tikhonov regularization matrix
   normalizeReg::AbstractRegularizationNormalization
 end
 
@@ -36,7 +34,6 @@ Creates a Kaczmarz object for the forward operator `A`.
 # Optional Keyword Arguments
   * `reg::AbstractParameterizedRegularization`          - regularization term
   * `normalizeReg::AbstractRegularizationNormalization` - regularization normalization scheme; options are `NoNormalization()`, `MeasurementBasedNormalization()`, `SystemMatrixBasedNormalization()`
-  * `weights::AbstractVector`                             - weights for the data term
   * `randomized::Bool`                                    - randomize Kacmarz algorithm
   * `subMatrixFraction::Real`                             - fraction of rows used in randomized Kaczmarz algorithm
   * `shuffleRows::Bool`                                   - randomize Kacmarz algorithm
@@ -46,25 +43,17 @@ Creates a Kaczmarz object for the forward operator `A`.
 See also [`createLinearSolver`](@ref), [`solve!`](@ref).
 """
 function Kaczmarz(A
-                ; reg = L2Regularization(0)
+                ; reg = L2Regularization(zero(real(eltype(A))))
                 , normalizeReg::AbstractRegularizationNormalization = NoNormalization()
-                , weights = nothing
                 , randomized::Bool = false
                 , subMatrixFraction::Real = 0.15
                 , shuffleRows::Bool = false
                 , seed::Int = 1234
                 , iterations::Int = 10
-                , regMatrix = nothing
                 )
 
   T = real(eltype(A))
 
-  # Apply Tikhonov regularization matrix
-  if regMatrix !== nothing
-    regMatrix = T.(regMatrix) # make sure regMatrix has the same element type as A
-    A = transpose(1 ./ sqrt.(regMatrix)) .* A # apply Tikhonov regularization to system matrix
-  end
-
   # Prepare regularization terms
   reg = isa(reg, AbstractVector) ? reg : [reg]
   reg = normalize(Kaczmarz, normalizeReg, reg, A, nothing)
@@ -76,6 +65,11 @@ function Kaczmarz(A
     deleteat!(reg, idx)
   end
 
+  # Tikhonov matrix is only valid with NoNormalization or SystemMatrixBasedNormalization
+  if λ(L2) isa Vector && !(normalizeReg isa NoNormalization || normalizeReg isa SystemMatrixBasedNormalization)
+    error("Tikhonov matrix for Kaczmarz is only valid with no or system matrix based normalization")
+  end
+
   indices = findsinks(AbstractProjectionRegularization, reg)
   other = AbstractRegularization[reg[i] for i in indices]
   deleteat!(reg, indices)
@@ -86,14 +80,13 @@ function Kaczmarz(A
   end
   other = identity.(other)
 
-
-  # make sure weights are not empty
-  w = (weights!=nothing ? weights : ones(T,size(A,1)))
-
   # setup denom and rowindex
-  denom, rowindex = initkaczmarz(A, λ(L2), w)
+  A, denom, rowindex = initkaczmarz(A, λ(L2))
   rowIndexCycle = collect(1:length(rowindex))
-  probabilities = T.(rowProbabilities(A, rowindex))
+  probabilities = eltype(denom)[]
+  if randomized
+    probabilities = T.(rowProbabilities(A, rowindex))
+  end
 
   M,N = size(A)
   subMatrixSize = round(Int, subMatrixFraction*M)
@@ -101,13 +94,13 @@ function Kaczmarz(A
   u  = zeros(eltype(A),M)
   x = zeros(eltype(A),N)
   vl = zeros(eltype(A),M)
-  εw = zeros(eltype(A),length(rowindex))
+  εw = zero(eltype(A))
   τl = zero(eltype(A))
   αl = zero(eltype(A))
 
   return Kaczmarz(A, u, L2, other, denom, rowindex, rowIndexCycle, x, vl, εw, τl, αl,
-                  T.(w), randomized, subMatrixSize, probabilities, shuffleRows,
-                  Int64(seed), iterations, regMatrix,
+                  randomized, subMatrixSize, probabilities, shuffleRows,
+                  Int64(seed), iterations,
                   normalizeReg)
 end
 
@@ -117,36 +110,46 @@ end
 (re-) initializes the Kacmarz iterator
 """
 function init!(solver::Kaczmarz, b; x0 = 0)
+  λ_prev = λ(solver.L2)
   solver.L2  = normalize(solver, solver.normalizeReg, solver.L2,  solver.A, b)
   solver.reg = normalize(solver, solver.normalizeReg, solver.reg, solver.A, b)
 
   λ_ = λ(solver.L2)
 
+  # λ changed => recompute denoms
+  if λ_ != λ_prev
+    # A must be unchanged, since we do not store the original SM
+    _, solver.denom, solver.rowindex = initkaczmarz(solver.A, λ_)
+    solver.rowIndexCycle = collect(1:length(rowindex))
+    if solver.randomized
+      solver.probabilities = T.(rowProbabilities(solver.A, rowindex))
+    end
+  end
+
   if solver.shuffleRows || solver.randomized
     Random.seed!(solver.seed)
   end
   if solver.shuffleRows
     shuffle!(solver.rowIndexCycle)
   end
-  solver.u .= b
 
   # start vector
   solver.x .= x0
   solver.vl .= 0
 
-  for i=1:length(solver.rowindex)
-    j = solver.rowindex[i]
-    solver.ɛw[i] = sqrt(λ_) / solver.weights[j]
+  solver.u .= b
+  if λ_ isa Vector
+    solver.ɛw = 0
+  else
+    solver.ɛw = sqrt(λ_)
   end
 end
 
-function solversolution(solver::Kaczmarz)
-  # backtransformation of solution with Tikhonov matrix
-  if solver.regMatrix !== nothing
-    return solver.x .* (1 ./ sqrt.(solver.regMatrix))
-  end
-  return solver.x
+
+function solversolution(solver::Kaczmarz{matT, RN}) where {matT, R<:L2Regularization{<:Vector}, RN <: Union{R, AbstractNestedRegularization{<:R}}}
+  return solver.x .* (1 ./ sqrt.(λ(solver.L2)))
 end
+solversolution(solver::Kaczmarz) = solver.x
 solverconvergence(solver::Kaczmarz) = (; :residual => norm(solver.vl))
 
 function iterate(solver::Kaczmarz, iteration::Int=0)
@@ -159,11 +162,8 @@ function iterate(solver::Kaczmarz, iteration::Int=0)
   end
 
   for i in usedIndices
-    j = solver.rowindex[i]
-    solver.τl = dot_with_matrix_row(solver.A,solver.x,j)
-    solver.αl = solver.denom[i]*(solver.u[j]-solver.τl-solver.ɛw[i]*solver.vl[j])
-    kaczmarz_update!(solver.A,solver.x,j,solver.αl)
-    solver.vl[j] += solver.αl*solver.ɛw[i]
+    row = solver.rowindex[i]
+    iterate_row_index(solver, solver.A, row, i)
   end
 
   for r in solver.reg
@@ -173,6 +173,14 @@ function iterate(solver::Kaczmarz, iteration::Int=0)
   return solver.vl, iteration+1
 end
 
+iterate_row_index(solver::Kaczmarz, A::AbstractLinearSolver, row, index) = iterate_row_index(solver, Matrix(A[row, :]), row, index) 
+function iterate_row_index(solver::Kaczmarz, A, row, index)
+  solver.τl = dot_with_matrix_row(A,solver.x,row)
+  solver.αl = solver.denom[index]*(solver.u[row]-solver.τl-solver.ɛw*solver.vl[row])
+  kaczmarz_update!(A,solver.x,row,solver.αl)
+  solver.vl[row] += solver.αl*solver.ɛw
+end
+
 @inline done(solver::Kaczmarz,iteration::Int) = iteration>=solver.iterations
 
 
@@ -180,41 +188,47 @@ end
 This function calculates the probabilities of the rows of the system matrix
 """
 
-function rowProbabilities(A::AbstractMatrix, rowindex)
-  M,N = size(A)
-  normS = norm(A)
+function rowProbabilities(A, rowindex)
+  normA² = rownorm²(A, 1:size(A, 1))
   p = zeros(length(rowindex))
   for i=1:length(rowindex)
     j = rowindex[i]
-    p[i] = (norm(A[j,:]))^2 / (normS)^2
+    p[i] = rownorm²(A, j) / (normA²)
   end
-
   return p
 end
 
+
 ### initkaczmarz ###
 
 """
-    initkaczmarz(A::AbstractMatrix,λ,weights::Vector)
+    initkaczmarz(A::AbstractMatrix,λ)
 
 This function saves the denominators to compute αl in denom and the rowindices,
 which lead to an update of x in rowindex.
 """
-function initkaczmarz(A::AbstractMatrix,λ,weights::Vector)
-  T = typeof(real(A[1]))
+function initkaczmarz(A,λ)
+  T = real(eltype(A))
   denom = T[]
   rowindex = Int64[]
 
-  for i=1:size(A,1)
-    s² = rownorm²(A,i)*weights[i]^2
+  for i = 1:size(A, 1)
+    s² = rownorm²(A,i)
     if s²>0
-      push!(denom,weights[i]^2/(s²+λ))
+      push!(denom,1/(s²+λ))
       push!(rowindex,i)
     end
   end
-  denom, rowindex
+  return A, denom, rowindex
+end
+function initkaczmarz(A, λ::Vector)
+  λ = real(eltype(A)).(λ)
+  A = initikhonov(A, λ)
+  return initkaczmarz(A, 0)
 end
 
+initikhonov(A, λ) = transpose((1 ./ sqrt.(λ)) .* transpose(A)) # optimize structure for row access
+initikhonov(prod::ProdOp{Tc, WeightingOp{T}, matT}, λ) where {T, Tc<:Union{T, Complex{T}}, matT} = ProdOp(prod.A, initikhonov(prod.B, λ))
 ### kaczmarz_update! ###
 
 """
@@ -242,6 +256,11 @@ function kaczmarz_update!(B::Transpose{T,S}, x::Vector,
   end
 end
 
+function kaczmarz_update!(prod::ProdOp{Tc, WeightingOp{T}, matT}, x::Vector, k, beta) where {T, Tc<:Union{T, Complex{T}}, matT}
+  weight = prod.A.weights[k]
+  kaczmarz_update!(prod.B, x, k, weight*beta) # only for real weights
+end
+
 # kaczmarz_update! with manual simd optimization
 for (T,W, WS,shufflevectorMask,vσ) in [(Float32,:WF32,:WF32S,:shufflevectorMaskF32,:vσF32),(Float64,:WF64,:WF64S,:shufflevectorMaskF64,:vσF64)]
     eval(quote
diff --git a/src/Regularization/NormalizedRegularization.jl b/src/Regularization/NormalizedRegularization.jl
index 97a40ebd..b8479a1b 100644
--- a/src/Regularization/NormalizedRegularization.jl
+++ b/src/Regularization/NormalizedRegularization.jl
@@ -30,6 +30,7 @@ Nested regularization term that scales `λ` according to normalization scheme. T
 struct NormalizedRegularization{T, S, R} <: AbstractScaledRegularization{T, S}
   reg::R
   factor::T
+  NormalizedRegularization(reg::R, factor) where {T, R <: AbstractParameterizedRegularization{<:AbstractArray{T}}} = new{T, R, R}(reg, factor)
   NormalizedRegularization(reg::R, factor) where {T, R <: AbstractParameterizedRegularization{T}} = new{T, R, R}(reg, factor)
   NormalizedRegularization(reg::R, factor) where {T, RN <: AbstractParameterizedRegularization{T}, R<:AbstractNestedRegularization{RN}} = new{T, RN, R}(reg, factor)
 end
@@ -43,17 +44,16 @@ normalize(::MeasurementBasedNormalization, A, b::Nothing) = one(real(eltype(A)))
 
 normalize(::SystemMatrixBasedNormalization, ::Nothing, _) = error("SystemMatrixBasedNormalization requires supplying A to the constructor of the solver")
 
-function normalize(::SystemMatrixBasedNormalization, A::AbstractArray{T}, b) where {T}
+function normalize(::SystemMatrixBasedNormalization, A, b)
   M = size(A, 1)
   N = size(A, 2)
 
-  energy = zeros(T, M)
+  energy = zeros(real(eltype(A)), M)
   for m=1:M
     energy[m] = sqrt(rownorm²(A,m))
   end
 
   trace = norm(energy)^2/N
-  # TODO where setlamda? here we dont know λ
   return trace
 end
 normalize(::NoNormalization, A, b) = nothing
diff --git a/src/Regularization/ScaledRegularization.jl b/src/Regularization/ScaledRegularization.jl
index eca47e8f..22f9fcce 100644
--- a/src/Regularization/ScaledRegularization.jl
+++ b/src/Regularization/ScaledRegularization.jl
@@ -6,7 +6,7 @@ Nested regularization term that applies a `scalefactor` to the regularization pa
 
 See also [`scalefactor`](@ref), [`λ`](@ref), [`innerreg`](@ref).
 """
-abstract type AbstractScaledRegularization{T, S<:AbstractParameterizedRegularization{T}} <: AbstractNestedRegularization{S} end
+abstract type AbstractScaledRegularization{T, S<:AbstractParameterizedRegularization{<:Union{T, <:AbstractArray{T}}}} <: AbstractNestedRegularization{S} end
 """
     scalescalefactor(reg::AbstractScaledRegularization)
 
@@ -20,7 +20,7 @@ return `λ` of `inner` regularization term scaled by `scalefactor(reg)`.
 
 See also [`scalefactor`](@ref), [`innerreg`](@ref).
 """
-λ(reg::AbstractScaledRegularization) = λ(innerreg(reg)) * scalefactor(reg)
+λ(reg::AbstractScaledRegularization) = λ(innerreg(reg)) .* scalefactor(reg)
 
 export FixedScaledRegularization
 struct FixedScaledRegularization{T, S, R} <: AbstractScaledRegularization{T, S}
diff --git a/src/Utils.jl b/src/Utils.jl
index 612d7d99..8b4d29f7 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -3,15 +3,15 @@ export rownorm², nrmsd
 """
 This function computes the 2-norm² of a rows of S for dense matrices.
 """
-function rownorm²(B::Transpose{T,S},row::Int) where {T,S<:DenseMatrix}
+function rownorm²(B::Transpose{T,S},row::Int64) where {T,S<:DenseMatrix}
   A = B.parent
-  U = typeof(real(A[1]))
+  U = real(eltype(A))
   res::U = BLAS.nrm2(size(A,1), pointer(A,(LinearIndices(size(A)))[1,row]), 1)^2
   return res
 end
 
-function rownorm²(A::AbstractMatrix,row::Int)
-  T = typeof(real(A[1]))
+function rownorm²(A::AbstractMatrix,row::Int64)
+  T = real(eltype(A))
   res = zero(T)
   @simd for n=1:size(A,2)
     res += abs2(A[row,n])
@@ -19,17 +19,26 @@ function rownorm²(A::AbstractMatrix,row::Int)
   return res
 end
 
+rownorm²(A::AbstractLinearOperator,row::Int64) = rownorm²(Matrix(A[row, :]), 1)
+rownorm²(A::ProdOp{T, <:WeightingOp, matT}, row::Int64) where {T, matT} = A.A.weights[row]^2*rownorm²(A.B, row)
+
 """
-This function computes the 2-norm² of a rows of S for dense matrices.
+This function computes the 2-norm² of a rows of S for sparse matrices.
 """
-function rownorm²(B::Transpose{T,S},row::Int) where {T,S<:SparseMatrixCSC}
+function rownorm²(B::Transpose{T,S},row::Int64) where {T,S<:SparseMatrixCSC}
   A = B.parent
-  U = typeof(real(A[1]))
+  U = real(eltype(A))
   res::U = BLAS.nrm2(A.colptr[row+1]-A.colptr[row], pointer(A.nzval,A.colptr[row]), 1)^2
   return res
 end
 
-
+function rownorm²(A, rows)
+  res = zero(real(eltype(A)))
+  @simd for row in rows
+    res += rownorm²(A, row)
+  end
+  return res
+end
 
 
 ### dot_with_matrix_row ###
@@ -90,6 +99,10 @@ function dot_with_matrix_row(B::Transpose{T,S},
   tmp
 end
 
+function dot_with_matrix_row(prod::ProdOp{T, <:WeightingOp, matT}, x::Vector{T}, k) where {T, matT}
+  A = prod.B
+  return prod.A.weights[k]*dot_with_matrix_row(A, x, k)
+end
 
 
 
diff --git a/src/proximalMaps/ProxL2.jl b/src/proximalMaps/ProxL2.jl
index 8e000354..fd1abe08 100644
--- a/src/proximalMaps/ProxL2.jl
+++ b/src/proximalMaps/ProxL2.jl
@@ -16,7 +16,7 @@ end
 performs the proximal map for Tikhonov regularization.
 """
 function prox!(::L2Regularization, x::AbstractArray{Tc}, λ::T) where {T, Tc <: Union{T, Complex{T}}}
-  x[:] .*= 1. / (1. + 2. *λ)#*x
+  x[:] .*= 1. ./ (1. .+ 2. .*λ)#*x
   return x
 end
 
@@ -26,3 +26,10 @@ end
 returns the value of the L2-regularization term
 """
 norm(::L2Regularization, x::AbstractArray{Tc}, λ::T) where {T, Tc <: Union{T, Complex{T}}} = λ*norm(x,2)^2
+function norm(::L2Regularization, x::AbstractArray{Tc}, λ::AbstractArray{T}) where {T, Tc <: Union{T, Complex{T}}}
+  res = zero(real(eltype(x)))
+  for i in eachindex(x)
+    res+= λ[i]*abs2(x[i])
+  end
+  return res
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index d1790b52..864b30b5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-using RegularizedLeastSquares, LinearAlgebra
+using RegularizedLeastSquares, LinearAlgebra, RegularizedLeastSquares.LinearOperatorCollection
 # Packages for testing only
 using Random, Test
 using FFTW
diff --git a/test/testKaczmarz.jl b/test/testKaczmarz.jl
index 880d1830..0e85c8c9 100644
--- a/test/testKaczmarz.jl
+++ b/test/testKaczmarz.jl
@@ -42,7 +42,7 @@ end
   regMatrix = rand(2) # Tikhonov matrix
 
   solver = Kaczmarz
-  S = createLinearSolver(solver, A, iterations=200, regMatrix=regMatrix)
+  S = createLinearSolver(solver, A, iterations=200, reg=[L2Regularization(regMatrix)])
   x_approx = solve!(S,b)
   #@info "Testing solver $solver ...: $x  == $x_approx"
   @test norm(x - x_approx) / norm(x) ≈ 0 atol=0.1
@@ -61,7 +61,7 @@ end
   # @show A, x, regMatrix
   # use regularization matrix
 
-  S = createLinearSolver(solver, A, iterations=100, regMatrix=regMatrix)
+  S = createLinearSolver(solver, A, iterations=100, reg=[L2Regularization(regMatrix)])
   x_matrix = solve!(S,b)
 
   # use standard reconstruction
@@ -73,6 +73,26 @@ end
   @test norm(x_approx - x_matrix) / norm(x_approx) ≈ 0 atol=0.1
 end
 
+@testset "Kaczmarz Weighting Matrix" begin
+  M = 12
+  N = 8
+  A = rand(M,N)+im*rand(M,N)
+  x = rand(N)+im*rand(N)
+  b = A*x
+  w = WeightingOp(rand(M))
+  d = diagm(w.weights)
+
+  reg = L2Regularization(rand())
+
+  solver = Kaczmarz
+  S = createLinearSolver(solver, d*A, iterations=200, reg = reg)
+  S_weighted = createLinearSolver(solver, *(ProdOp, w, A), iterations=200, reg = reg)
+  x_approx = solve!(S, d*b)
+  x_weighted = solve!(S_weighted, d*b)
+  #@info "Testing solver $solver ...: $x  == $x_approx"
+  @test isapprox(x_approx, x_weighted)
+end
+
 
 # Test Kaczmarz parameters
 @testset "Kaczmarz parameters" begin