JuliaSmoothOptimizers / Krylov.jl

A Julia Basket of Hand-Picked Krylov Methods
Other
349 stars 51 forks source link

[documentation] Sensitivity analysis #839

Open amontoison opened 10 months ago

amontoison commented 10 months ago

I should add an example in the documentation to explain how we can do sensitivity analysis with Krylov.jl.

using ChainRulesCore
using Krylov
using Test
using LinearAlgebra
using SparseArrays
using ForwardDiff
import ForwardDiff: Dual, partials, value

include("test/test_utils.jl")

A, b = sparse_laplacian(FC=Float64)

function Krylov.cg(_A::SparseMatrixCSC{Dual{T, V, N}, Int64}, _b::Vector{Dual{T, V, N}}; options...) where {T, V, N}
  A = SparseMatrixCSC(_A.m, _A.n, _A.colptr, _A.rowval, value.(_A.nzval))
  dA = SparseMatrixCSC(_A.m, _A.n, _A.colptr, _A.rowval, partials.(_A.nzval, 1))
  b = value.(_b)
  db = partials.(_b,1)
  x, stats = cg(A,b)
  nb = db - dA*x
  dx, dstats = cg(dA,nb)
  return (Dual.(x, dx), stats)
end

x, stats = cg(A,b)

dA = ForwardDiff.Dual.(A, A)
db = ForwardDiff.Dual.(b, b)

It solves the normal value system Ax = b and the tangent system A*dx + dA*x= db (or A*dx = b - dA*x).