Efficient Kaczmarz
Unlike many of the other solvers provided by RegularizedLeastSquares.jl, the Kaczmarz method does not utilize a matrix-vector product with the operator $\mathbf{A}$ nor the normal operator $\mathbf{A*A}$. Instead, it uses the rows of $\mathbf{A}$ to update the solution iteratively. Efficient Kaczmarz implementation therefore require very efficient dot products with the rows of $\mathbf{A}$. In RegularizedLeastSquares.jl, this is achieved with the dot_with_matrix_row function.
using RegularizedLeastSquares
A = randn(256, 256)
x = randn(256)
b = A*x;The dot_with_matrix_row function calculates the dot product between a row of A and the current approximate solution of x:
row = 1
isapprox(RegularizedLeastSquares.dot_with_matrix_row(A, x, row), sum(A[row, :] .* x))trueSince in Julia, dense arrays are stored in column-major order, such a row-based operation is quite inefficient. A workaround is to transpose the matrix then pass it to a Kaczmarz solver.
At = collect(transpose(A))
A_eff = transpose(At)256×256 transpose(::Matrix{Float64}) with eltype Float64:
-0.275864 -0.73746 0.853414 … 0.293825 -0.502286 -0.11887
0.898753 0.137198 0.998105 0.407286 -1.22362 0.610509
-0.721099 1.38746 -1.01862 0.879271 -2.34228 -0.523434
0.7809 -0.887607 -2.03852 2.23992 -0.187041 -2.4077
0.971881 2.20801 -0.0103514 0.172781 -0.37087 1.52988
-0.485217 0.425664 0.0914575 … 0.55081 0.493884 0.897434
-0.710359 0.851985 0.4346 -0.877772 1.74841 -0.176348
-2.01853 -0.987419 -0.5984 0.819381 -1.53312 0.853734
0.870522 -0.0517138 -0.29153 0.28371 -2.23615 0.718148
0.447491 -1.07847 1.49144 0.0516904 2.82065 0.682398
⋮ ⋱ ⋮
-1.52396 0.729011 -0.151012 -1.20962 0.807203 -1.35141
1.36549 -1.6457 1.87423 0.183826 -0.568891 -0.0148501
-0.728423 2.46315 1.11471 1.26203 -0.390659 1.13454
0.376685 1.19146 0.87661 … -0.361811 -1.12049 -1.98131
1.58625 1.25387 -0.256293 -2.79249 -0.871554 0.188822
0.692525 0.0836304 -0.730701 -0.99716 0.220903 2.52113
-0.492048 -0.301455 0.0955891 -0.375775 -0.118445 1.89548
1.29308 -0.107961 1.27551 -0.973497 2.13307 -0.215691
0.361892 -0.287885 0.218406 … 0.811435 -0.916367 2.45249Note that the transpose function can return a lazy transpose object, so we first collect the transpose into a dense matrix. Then we transpose it again to get the efficient representation of the matrix.
We can compare the performance using the BenchmarkTools.jl package. First for the original matrix:
using BenchmarkTools
solver = createLinearSolver(Kaczmarz, A; reg = L2Regularization(0.0001), iterations=100)
@benchmark solve!(solver, b) samples = 100BenchmarkTools.Trial: 100 samples with 1 evaluation per sample.
Range (min … max): 27.751 ms … 28.353 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 28.167 ms ┊ GC (median): 0.00%
Time (mean ± σ): 28.171 ms ± 54.119 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▃█
▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▆▄███▆▆▆▃▇▄▁▁▃▁▁▃ ▃
27.8 ms Histogram: frequency by time 28.3 ms <
Memory estimate: 17.30 KiB, allocs estimate: 504.And then for the efficient matrix:
solver_eff = createLinearSolver(Kaczmarz, A_eff; reg = L2Regularization(0.0001), iterations=100)
@benchmark solve!(solver_eff, b) samples = 100BenchmarkTools.Trial: 100 samples with 1 evaluation per sample.
Range (min … max): 1.402 ms … 1.585 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 1.417 ms ┊ GC (median): 0.00%
Time (mean ± σ): 1.423 ms ± 30.186 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▁▃▅█▅▅
▇██████▄▄▄▃▄▃▃▁▁▁▁▁▁▁▃▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▃ ▃
1.4 ms Histogram: frequency by time 1.58 ms <
Memory estimate: 17.31 KiB, allocs estimate: 505.We can also combine the efficient matrix with a weighting matrix, as is shown in the Weighting example.
Custom operators need to implement the dot_with_matrix_row function to be used with the Kaczmarz solver. Ideally, such an implementation is allocation free.
This page was generated using Literate.jl.