Jared McBride

Jared McBride

I am a PhD candatite in the University of Arizona Aplied Math GIDP

Comparing REK to QR

REK vs QR

REK v. QR

In [83]:
using JLD
using PyPlot
using LinearAlgebra, Statistics

rek = include("../../../../Kaczmarz/Kaczmarz.jl")
util = include("../../../Tools/Utils.jl")
mrl = include("../../../Tools/WFMR_LS.jl")
mrk = include("../../../Tools/WFMR_kac.jl")

at = include("../../../Tools/AnalysisToolbox.jl")

WARNING: replacing module Kaczmarz.
WARNING: replacing module VariousUtilities.
WARNING: replacing module WFMR_LS.
WARNING: replacing module WFMR_Kac.
WARNING: replacing module AnalysisToolbox.
Out[83]:
Main.AnalysisToolbox
In [84]:
# Load Data

#This is writen to allow for loading any of 5 different runs of the KSE solution


N = 499999


geni = 1
gen = "lin1e5_r$geni"     # this is just a reference designation it shows up in the
                          # output file. I think of generatrion.

savefile_pred = "../../../../data/KSE_Data/ks_pred_$gen.jld"
println("Sol load location: " * savefile_pred)

Pred = load(savefile_pred,"pred");
Sig  = load(savefile_pred,"sig")


sig = Sig[1:1,1:N]
pred = Pred[1:1,1:N]

Sol load location: ../../../../data/KSE_Data/ks_pred_lin1e5_r1.jld
Out[84]:
1×499999 Matrix{ComplexF64}:
 -0.0643694-0.0698753im  -0.0622791-0.0665026im  …  0.13453+0.11138im
In [85]:
M_out = 20
verb = true
vew = 1000:2000

# Form design matrix
O = @timed mrl.get_hash(pred; M_out)
Xh = O.value
verb && println("Design matrix formed. It took $(O.time) sec")

# Extract and reshape signal
Y = Array(transpose(sig[:,M_out:end]));

Design matrix formed. It took 0.173529269 sec
In [86]:
h_bs = Xh[1:5000,:] \ Y[1:5000]

Y_hat = Xh * h_bs

res = util.TakeLook(Y,Y_hat; vew)
println("MSE: ",var(res))

MSE: 3.0834285029278176e-23
In [87]:
h_rek = rek.solve(Xh[1:5000,:], Y[1:5000], eps = 1e-14, maxcount = 1000, reportperiod = 60, verbose = verb).sol

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

(m, n, subcount, Asum) = (5000, 20, 160, 4746.468002173153)
#Kaczmarz: 1000 outer loops reached
MSE: 1.026875486466541e-6

Preconditioning (Cholesky factor)

In [118]:
function precond(Xh,pcond;verb = true, eps = 1e-6)
    M = M_out
    w = I

    for i = 1:pcond
        σ = eps^2*var(Xh[:,M])
        A = Xh'Xh/size(Xh,1) + σ*I
        verb && println("preconditioning $i original condition number: ",cond(A))
        L = cholesky(A).L
        verb && println("preconditioning $i Cholesky condition number: ",cond(L))

        W = inv(conj(L))
        w = W*w

        Xh = Xh*transpose(W)
        verb && println("preconditioning $i new condition number: ",cond(Xh))
    end
    Xh, w
end
Out[118]:
precond (generic function with 1 method)
In [138]:
Xhw, w = precond(Xh, 1)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-14, maxcount = 1000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
(m, n, subcount, Asum) = (5000, 20, 160, 43765.41247400333)
#Kaczmarz: 1000 outer loops reached
outercount = 1000
innercount = 160000
MSE: 7.481775074448454e-11
In [136]:
Xhw, w = precond(Xh, 2)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-14, maxcount = 10000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
preconditioning 2 original condition number: 4.353725483168935e12
preconditioning 2 Cholesky condition number: 2.086528288515005e6
preconditioning 2 new condition number: 988.5476925920568
(m, n, subcount, Asum) = (5000, 20, 160, 91501.8725819039)
#Kaczmarz: 10000 outer loops reached
outercount = 10000
innercount = 1600000
MSE: 1.3308198589604643e-11
In [139]:
Xhw, w = precond(Xh, 3)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-14, maxcount = 10000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
preconditioning 2 original condition number: 4.353725483168935e12
preconditioning 2 Cholesky condition number: 2.086528288515005e6
preconditioning 2 new condition number: 988.5476925920568
preconditioning 3 original condition number: 977226.0537615535
preconditioning 3 Cholesky condition number: 988.5474463705507
preconditioning 3 new condition number: 1.0000002490739892
(m, n, subcount, Asum) = (5000, 20, 160, 117413.78115309498)
#Kaczmarz: early exit
outercount = 8
innercount = 1280
MSE: 1.4746741716675059e-18
In [141]:
Xhw, w = precond(Xh, 3)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-16, maxcount = 10000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
preconditioning 2 original condition number: 4.353725483168935e12
preconditioning 2 Cholesky condition number: 2.086528288515005e6
preconditioning 2 new condition number: 988.5476925920568
preconditioning 3 original condition number: 977226.0537615535
preconditioning 3 Cholesky condition number: 988.5474463705507
preconditioning 3 new condition number: 1.0000002490739892
(m, n, subcount, Asum) = (5000, 20, 160, 117413.78115309498)
#Kaczmarz: early exit
outercount = 247
innercount = 39520
MSE: 4.410457448010638e-18
In [145]:
Xhw, w = precond(Xh, 3)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-18, maxcount = 20000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
preconditioning 2 original condition number: 4.353725483168935e12
preconditioning 2 Cholesky condition number: 2.086528288515005e6
preconditioning 2 new condition number: 988.5476925920568
preconditioning 3 original condition number: 977226.0537615535
preconditioning 3 Cholesky condition number: 988.5474463705507
preconditioning 3 new condition number: 1.0000002490739892
(m, n, subcount, Asum) = (5000, 20, 160, 117413.78115309498)
#Kaczmarz: 20000 outer loops reached
outercount = 20000
innercount = 3200000
MSE: 8.194925747719924e-18
In [147]:
Xhw, w = precond(Xh, 3)

h_rek, outercount, innercount = rek.solve(Xhw[1:5000,:], Y[1:5000],
                                    eps = 1e-20, maxcount = 200000, reportperiod = 60, verbose = verb)
@show outercount
@show innercount

h_rek = transpose(w)*h_rek

Y_hat_rek = Xh * h_rek

res = util.TakeLook(Y,Y_hat_rek; vew)
println("MSE: ",var(res))

preconditioning 1 original condition number: 1.9799050040923062e13
preconditioning 1 Cholesky condition number: 4.449514319348537e6
preconditioning 1 new condition number: 2.0610299857304611e9
preconditioning 2 original condition number: 4.353725483168935e12
preconditioning 2 Cholesky condition number: 2.086528288515005e6
preconditioning 2 new condition number: 988.5476925920568
preconditioning 3 original condition number: 977226.0537615535
preconditioning 3 Cholesky condition number: 988.5474463705507
preconditioning 3 new condition number: 1.0000002490739892
(m, n, subcount, Asum) = (5000, 20, 160, 117413.78115309498)
#Kaczmarz: 7510932/32000000 steps took 1 min 0.0 sec; eta 3 min 15.627 sec
#Kaczmarz: 15004481/32000000 steps took 2 min 0.0 sec; eta 2 min 15.923 sec
#Kaczmarz: 22496977/32000000 steps took 3 min 0.0 sec; eta 1 min 16.034 sec
#Kaczmarz: 29997423/32000000 steps took 4 min 0.0 sec; eta 16.021 sec

#Kaczmarz: 200000 outer loops reached
outercount = 200000
innercount = 32000000
MSE: 6.086692244477209e-19
In [ ]: