Sullivan p. 96

bayesian_regularization.jl

let
	function ℓπ(u,A,y,λ₁=1e-4,λ₂=1e-4) 
	        n,m = size(A)
	        AQiA =  A'*(I./λ₁)*A
	        P = AQiA+I./λ₁
	        K = pinv(AQiA)*A'*I./λ₁
	        [-1.0.*(-0.5*sum((u.-inv(P)*(AQiA*K*y+ I./λ₂*u)).^2))]
	    end
		rng = Xoshiro(123)
		A = rand(rng,5,5)
		x = rand(rng,5,1)
	
		y = A*x 
		u0 = rand(rng,5)
		f_opti(x,p) = ℓπ(x,A,y)
	    prob = NonlinearProblem(NonlinearFunction(f_opti), u0, 0.0)
	    sol = solve(prob, NonlinearSolve.LevenbergMarquardt())
		@info sol
		x_new = sol.u[:,:]
		@info x
		@info x_new
		@info sum(A*x.-y)
		@info sum(A*x_new.-y)
end