GT ICR - Nicolas Monnier

Transcript

1. L’imagerie en radio interférométrie : Réduire le coût calculatoire avec la méthode Grid to Grid Présentation pour le GT-ICR Nicolas Monnier 28 Avril 2025
2. Table of contents 1 Radio Interferometric Imaging 2 Grid to Grid (G2G) Nicolas Monnier 2
3. Table of contents 1 Radio Interferometric Imaging 2 Grid to Grid (G2G) Nicolas Monnier 3
4. Radio interferometry Radio Interferometer : Array of antennas. Angular resolution : ~ Nicolas Monnier 4
5. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
6. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] Visibility (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
7. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] True Sky (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
8. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] Direction Independent Effects (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
9. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] Direction Dependent Effects (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
10. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] Observation Delay (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
11. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Nicolas Monnier 5
12. Radio Interferometric Measurement Equation (RIME) Jones formulation [Hamaker, 1996; Smirnov, 2011] (l, m, n): Direction of observation coordinate system (u, v, w): Visibility coordinate system Simplified formulation vpq(u, v) = lm x(l, m)e−2iπ(ul+vm) dl dm Nicolas Monnier 5
13. Discretization of the RIME Forward problem SKA-Mid SKA-Low 197 Dishes 110.000 Antennas v ∈ CM : Visibility vector x ∈ RP : Sky vector n ∈ CM : Gaussian noise Nicolas Monnier 6
14. Discretization of the RIME Forward problem SKA-Mid SKA-Low 197 Dishes 110.000 Antennas v ∈ CM : Visibility vector x ∈ RP : Sky vector n ∈ CM : Gaussian noise Instrumental Model H • Irregular sampling of the Fourier space : NUFFT • Prohibitive computational cost • Approximated by Fast Fourier Transform and sampling Nicolas Monnier 6
15. Discretization of the RIME Forward problem v ∈ CM : Visibility vector x ∈ RP : Sky vector n ∈ CM : Gaussian noise F ∈ CP ×P : Fourier transform matrix S ∈ CM×P : Sampling operator (Degridding) Instrumental Model H • Irregular sampling of the Fourier space : NUFFT • Prohibitive computational cost • Approximated by Fast Fourier Transform and sampling Nicolas Monnier 6
16. Synthesis imaging y ∈ CP : Dirty image F † ∈ CP ×P : Fourier transform matrix S† ∈ CP ×M : Sampling operator (Gridding) v ∈ CM : Visibility vector Nicolas Monnier 7
17. uv-coverage VLA antennas position uv-coverage - 10min Dirty image Nicolas Monnier 8
18. uv-coverage VLA antennas position uv-coverage - 4h Dirty image Nicolas Monnier 8
19. uv-coverage VLA antennas position uv-coverage - 8h Dirty image Nicolas Monnier 8
20. Fourier synthesis problem • Incomplete sampling of the Fourier space v = SF x + n • Ill-posed inverse problem : infinity of solutions • Solution defines as minimizing a mix-criterion : well-posed problem x = arg min x ∥v − SF x∥2 2 + αR(x) • Iterative algorithms involve computation of the data fidelity gradient ∇J(x) = F †S†(v − SF x) Nicolas Monnier 9
21. Iterative algorithm - Nicolas Monnier 10
22. Iterative algorithm - Compute gradient (Bottleneck) Nicolas Monnier 10
23. Iterative algorithm - Image Domain Gridder Nifty Gridder W-projection ... [Van der Tol, 2018] [Ye, 2022] [Cornwell, 2008] Compute gradient (Bottleneck) Nicolas Monnier 10
24. Iterative algorithm - Generate next estimate Image Domain Gridder Nifty Gridder W-projection ... [Van der Tol, 2018] [Ye, 2022] [Cornwell, 2008] Compute gradient (Bottleneck) Nicolas Monnier 10
25. Iterative algorithm - Generate next estimate Minor Cyles CLEAN Image Domain Gridder Nifty Gridder W-projection ... [Van der Tol, 2018] [Ye, 2022] [Cornwell, 2008] [Rau, 2011] [Högbom, 1974] [Schwab , 1984] Compute gradient (Bottleneck) Nicolas Monnier 10
26. Iterative algorithm - Generate next estimate CASA DDFacet WSClean Purify .... Minor Cyles CLEAN Image Domain Gridder Nifty Gridder W-projection ... [Van der Tol, 2018] [Ye, 2022] [Cornwell, 2008] [Rau, 2011] [Högbom, 1974] [Schwab , 1984] [Jaeger, 2008] [Tasse, 2017] [Offringa, 2014] [Carillo, 2014] Compute gradient (Bottleneck) Nicolas Monnier 10
27. Table of contents 1 Radio Interferometric Imaging 2 Grid to Grid (G2G) Nicolas Monnier 11
28. Gridding S† and degridding S operators Gridding Degridding Nicolas Monnier 12
29. Gradient computation optimization • High computing cost for Gridding and Degridding operators O(S†) = O(S) = MC2 supp • Gradient computation : ∇JST D (x) = F †S†(v − SF x) (1) ∇JG2G (x) = F †S†v − F †S†SF x (2) Fuse of S†S • Fast Holographic Deconvolution : H = S†S ∈ CP ×P , [Sullivan, 2012] • G2G : Reduce computational cost • G2G : Reduce memory footprint Nicolas Monnier 13
30. Gridding operator decomposition S† = A†C†O† • Accumulation A† on Oversampled Fourier grid • Discrete convolution C† of KNp × KNp grid • Sub-sampling O† Nicolas Monnier 14
31. Degridding operator decomposition S = ACO • Oversampling operator O • Discrete convolution C† of KNp × KNp grid • Mapping operator A Nicolas Monnier 15
32. Grid to Grid method (G2G) [Monnier et al., 2022] Fusion of sub-operators A†A • A∗ is diagonal matrix • M′ non-zero coefficient in A∗ • M′ ≪ M • ∇JG2G(x) = ∇JST D(x) Nicolas Monnier 16
33. G2G details • Implicit algorithm : no oversampled grid • Optimized GPU and multi-core CPU implementation • Precomputed convolution kernels • Method well-suited for compact configuration array • Generalisation of Baseline-Dependent Averaging (BDA) [Wijnholds, 2018; Atemkeng, 2016] Nicolas Monnier 17
34. Experiments VLA test case Dataset : 7.7 GB - Image size : 1280 × 1280 pixels Mean Absolute Error MeanAE(K) = 1 P P i=0 |∇J(xi,K=63 ) − ∇J(xi,K )| Nicolas Monnier 18
35. G2G implementations GPU implementation • GPU implementation based on [Romein, 2012] [Merry, 2016] • Minimize device memory access • Shared memory reduction CPU and GPU comparison • CPU : Dual-socket Ice Lake, 64 cores • GPU : NVIDIA A100-SXM • GPU 10× faster than CPU Nicolas Monnier 19
36. Memory footprint Compression ratio cα = M′ M Low cα is lower memory footprint and lower computational cost to compute ∇G2GJ(x) GB GB GB Nicolas Monnier 20
37. STD vs G2G comparison G2G STD/IDG Nicolas Monnier 21
38. GPU speedup State-of-the-art speedup comparison • Correction for non-coplanar array • STD : Standard w-projection gridder + degridder [Cornwell, 2008; Romein, 2012] • IDG : Image Domain Gridding [Van der Tol, 2018] Speedup = Time STD/IDG Time G2G 10 20 30 40 50 60 Oversampling factor K 0 25 50 75 100 125 Speedup for Csupp = 8 G2G/IDG G2G/STD Computation time : 10−2 –100 sec 10 20 30 40 50 60 Oversampling factor K 10 20 30 40 Speedup for Csupp = 64 G2G/IDG G2G/STD Computation time : 100 –101 sec Nicolas Monnier 22
39. GPU implementation comparison Comparison metric • Fourier-point throughput to compute g from g • MFT : Million Fourier-point Throughput (in FP/s) MFT = #Fourier point Time (S† + S) × 10−6 Nicolas Monnier 23
40. G2G Memory footprint K 8 16 24 32 40 64. STD single 7.7 7.7 7.7 7.7 7.7 7.7 double 15.4 15.4 15.4 15.4 15.4 15.4 G2G single 1.7 1.85 1.99 2.13 2.27 2.66 64 slices double 3.27 3.42 3.56 3.70 3.83 4.22 G2G single 0.17 0.32 0.46 0.60 0.73 1.12 1 slice double 0.19 0.34 0.48 0.62 0.76 1.15 Memory footprint needed to compute the gradient in GB. Nicolas Monnier 24