oclib: out of core optimization

1. oclib: out of core optimization Paul Sava paul@sep.stanford.edu

2. My goal • Quick introduction to the library • Not much about inversion theory • Ask the right questions paul@sep.stanford.edu

3. Jon Claerbout’s optimization book • Image estimation by example: Geophysical sounding image construction Multidimensional autoregression http://sepwww.stanford.edu/sep/prof paul@sep.stanford.edu

4. Linear operators L: linear operator m: image d: data paul@sep.stanford.edu

5. Linear operators m d L call operator_init(…) stat=operator(adj,add,m,d) paul@sep.stanford.edu

6. Least-squares inversion paul@sep.stanford.edu

7. Simple inversion m: model d: data L: physics R=Lm0-d iterate { g = L* R G = L g (m,R) <- step(m,R,g,G) } paul@sep.stanford.edu

8. Regularized inversion W: weighting A: regularization e : relative weight paul@sep.stanford.edu

9. Preconditioned inversion p: preconditioned m paul@sep.stanford.edu

10. Why out-of-core optimization Least-squares imaging m: image d: data L: WE modeling operator A: styling operator W: weighting operator Velocity analysis m: slowness perturbation d: image perturbation L: WE MVA operator A: styling operator W: weighting operator paul@sep.stanford.edu

11. Linear operator example: WE MVA call weimva_init( SLin=weiop_sllN_init, SCin=weiop_bor1_init, WCin=weiop_mwcN_init, FKin=weiop_wem3_init, FXin=weiop_ssf3_init, IGin=weiop_hcig_init) stat = weimva(adj,add, dS,dR, SLop=weiop_sllN, SCop=weiop_bor1, WCop=weiop_mwcN, FKop=weiop_wem3, FXop=weiop_ssf3, IGop=weiop_hcig) paul@sep.stanford.edu

12. Outline • Overview • Solvers • Examples paul@sep.stanford.edu

13. Simple solver init call oc_mpi_solver_init ( impi, ocniter, ocmaxmem, ocverb, mmovie, dmovie ) R=Lm-d iterate { g=L* R G=L g (m,R) <- step(m,R,g,G) } paul@sep.stanford.edu

14. Simple solver call oc_mpi_solver ( L=weimva, m_=dS, d_=dR, S=oc_cgstep, op1=weiop_slo1, op2=weiop_mwc1, op3=weiop_bor1, op4=weiop_wem3, op5=weiop_ssf3, op6=weiop_hcig ) R=Lm-d iterate { g=L* R G=L g (m,R) <- step(m,R,g,G) } paul@sep.stanford.edu

15. Simple solver R=Lm-d iterate { g=L* R G=L g (m,R) <- step(m,R,g,G) } do i=1,ocniter stat=L( T,F, g_, rd_, op1 … op9) stat=L( F,F, g_, gd_, op1 … op9) stat=S(forget,m_,g_,(/rd_/),(/gd_/),…) end do paul@sep.stanford.edu

16. Cluster execution R=Lm-d iterate { g=L* R G=L g (m,R) <- step(m,R,g,G) } do i=1,ocniter stat=L( T,F, g_, rd_, op1 … op9) stat=L( F,F, g_, gd_, op1 … op9) stat=S(forget,m_,g_,(/rd_/),(/gd_/),…) end do node 1 node 4 node 2 node 3 master paul@sep.stanford.edu

17. Solver with regularization call oc_mpi_solverreg( L=weimva, m_=dS, d_=dR, A=oc_laplacian, S=oc_cgstep, op1=weiop_slo1, op2=weiop_mwc1, op3=weiop_bor1, op4=weiop_wem3, op5=weiop_ssf3, op6=weiop_hcig ) paul@sep.stanford.edu

18. Solver with preconditioning call oc_mpi_solverpre( L=weimva, m_=dS, d_=dR, P=oc_ilaplacian, S=oc_cgstep, op1=weiop_slo1, op2=weiop_mwc1, op3=weiop_bor1, op4=weiop_wem3, op5=weiop_ssf3, op6=weiop_hcig ) paul@sep.stanford.edu

19. Outline • Overview • Solvers • Examples paul@sep.stanford.edu

20. Example & Lab • WE MVA for diffracted data • d: image perturbation • m: slowness perturbation • Operator calls • forward: scattering mode • adjoint: backprojection mode • Regularized inversion • model movie • data residual movie paul@sep.stanford.edu

21. Example paul@sep.stanford.edu

22. Example: d = Lm paul@sep.stanford.edu

23. Example: m = L*d paul@sep.stanford.edu

24. Example: m = [L*L+e2A*A]-1L*d paul@sep.stanford.edu

25. Iteration 1:Rd=d-Lm paul@sep.stanford.edu

26. Iteration 2:Rd=d-Lm paul@sep.stanford.edu

27. Iteration 3:Rd=d-Lm paul@sep.stanford.edu

28. Iteration 4:Rd=d-Lm paul@sep.stanford.edu

29. Iteration 5:Rd=d-Lm paul@sep.stanford.edu

30. Summary • Flexible, reusable f90 code • Cluster ready • Standard operator interface • Conjugate-gradient solvers • Regularization • Preconditioning paul@sep.stanford.edu

31. Resources (name@sep.stanford.edu) • Developers: • Paul, Bob • Manual • Marie • Users • Antoine, Nick paul@sep.stanford.edu