Extreme scale matrix factorizations in Exploration Seismology Felix - PowerPoint PPT Presentation

Extreme scale matrix factorizations in Exploration Seismology Felix J. Herrmann SLIM Georgia Institute of Technology Saturday, November 11, 17

SLIM Seismic inversion Infer 3D images & velocity models from mul+-experiment data: ‣ unknowns O (10 9 ) ‣ datapoints O (10 15 ) ‣ propagate wavelengths O (10 2 ) Saturday, November 11, 17

Wave-equation based inversion Retrieve medium parameters m from partial measurements of the solution of the wave-equation: A ( m ) u i = q i P i u i = d i q i m Saturday, November 11, 17

Wave-equation based inversion Large-scale parameter estimation problem: observed data M Φ ( m ) = 1 φ i ( m ) = 1 X 2 k P i u i � d i k 2 minimize M m i =1 ‣ number of field experiments is large O (10 3 − 10 5 ) M ‣ expensive to collect data points at total survey costs of $30 – 200 M O (10 6 − 10 7 ) d i ‣ expensive to evaluate flops per experiment w/ HPC costs of $25 – 500 M O (10 14 ) φ i ‣ is extremely large requiring local (= gradient-based ) optimization O (10 6 − 10 9 ) m Saturday, November 11, 17

Research questions “How can we exploit low-rank structure underlying surface seismic data and subsurface image volumes at low to mid frequencies?” ‣ reduce acquisition time, costs, environmental imprint of via randomized sampling & full azimuth processing ‣ lower storage & IO cost of wave-equation based inversion via on-the- fly data generation from data represented in factorized form ‣ form & manipulate massive full-subsurface offset image volumes via randomized probing of the double two-way wave equation Saturday, November 11, 17

Motivation – seismic surface data Large 5D volumes of seismic data 100’s of thousands of shots 6 Saturday, November 11, 17

Motivation – seismic surface data Large 5D volumes of seismic data 100’s of thousands of shots Will soon reach petabytes. 6 Saturday, November 11, 17

Motivation – image volumes Extremely large 6D image volumes quadratic in image size 7 Saturday, November 11, 17

Motivation – image volumes Extremely large 6D image volumes quadratic in image size Can not be formed explicitly on any computer... 7 Saturday, November 11, 17

Full-azimuth seismic data processing with coil acquisition Rajiv Kumar, Nick Moldoveanu, Keegan Lensink, and Felix J. Herrmann SLIM Georgia Institute of Technology Saturday, November 11, 17

3D seismic Challenge seismic data collected in 5 dimensions ‣ 1 for time ‣ 2 for the receivers ‣ 2 for the sources Compressive Sensing works well for vectorial transform-domain sparsity ‣ curvelets & other non-separable transforms are too slow & memory intensive ‣ prohibits scale up to 5D Can we exploit a different kind of structure ... Saturday, November 11, 17

Quick recap––matrix completion Aleksandr Y. Aravkin, Rajiv Kumar, Hassan Mansour, Ben Recht, and Felix J. Herrmann, “ Fast methods for denoising matrix completion formulations, with applications to robust seismic data interpolation ”, SIAM Journal on Scientific Computing , vol. 36, p. S237-S266, 2014 Rajiv Kumar, Haneet Wason, and Felix J. Herrmann, “ Source separation for simultaneous towed-streamer marine acquisition –- a compressed sensing approach ”, Geophysics , vol. 80, p. WD73-WD88, 2015. Rajiv Kumar, Curt Da Silva, Okan Akalin, Aleksandr Y. Aravkin, Hassan Mansour, Ben Recht, and Felix J. Herrmann, “ Efficient matrix completion for seismic data reconstruction ”, Geophysics , vol. 80, p. V97-V114, 2015. Curt Da Silva and Felix J. Herrmann, “ Optimization on the Hierarchical Tucker manifold - applications to tensor completion ”, Linear Algebra and its Applications , vol. 481, p. 131-173, 2015. 10 Saturday, November 11, 17

[Candes and Plan 2010, Oropeza and Sacchi 2011] Matrix completion ! ‣ signal structure ! - low rank/fast decay of singular values ! ! ‣ sampling scheme ! - missing data increase rank in “transform domain” ! ! ‣ recovery using rank penalization scheme Saturday, November 11, 17

Low-rank structure conventional 5D data, 5 Hz monochromatic slice, Sx-Sy matricization 10 20 30 40 50 60 100 70 80 200 90 100 Rec x, Rec y 300 0 10 Full data − 1 10 400 − 2 10 Normalized singular value − 3 10 500 − 4 10 600 − 5 10 100 200 300 400 500 600 − 6 10 Src x, Src y − 7 10 0 100 200 300 400 500 600 700 12 Saturday, November 11, 17

Low-rank structure conventional 5D data, 5 Hz monochromatic slice, Sx-Rx matricization 10 20 30 40 50 60 100 70 80 200 90 100 Src y, Rec y 300 0 10 (Rx Ry) matricization (Sy Ry) matricization − 1 10 400 − 2 10 Normalized singular value − 3 10 500 − 4 10 600 − 5 10 100 200 300 400 500 600 − 6 10 Src x, Rec x − 7 10 0 100 200 300 400 500 600 700 13 Saturday, November 11, 17

Matrix completion ! ‣ signal structure ! - low rank/fast decay of singular values ! ! ‣ sampling scheme ! - missing data increase rank in “transform domain” ! ! ‣ recovery using rank penalization scheme Saturday, November 11, 17

Low-rank structure jittered data, 5 Hz monochromatic slice, Sx-Sy matricization 10 20 30 40 50 60 100 70 80 200 90 100 Rec x, Rec y 300 0 10 No subsampling 50% missing sources − 1 400 10 − 2 10 Normalized singular value 500 − 3 10 − 4 10 600 − 5 10 100 200 300 400 500 600 Src x, Src y − 6 10 − 7 10 0 100 200 300 400 500 600 700 15 Saturday, November 11, 17

Low-rank structure jittered data, 5 Hz monochromatic slice, Sx-Rx matricization 10 20 30 40 50 60 100 70 80 200 90 100 Src y, Rec y 300 0 10 No subsampling 50% missing sources − 1 10 400 − 2 10 Normalized singular value − 3 10 500 − 4 10 − 5 10 600 − 6 100 200 300 400 500 600 10 Src x, Rec x − 7 10 − 8 10 0 100 200 300 400 500 600 700 16 Saturday, November 11, 17

Matrix completion ! ‣ signal structure ! - low rank/fast decay of singular values ! ! ‣ sampling scheme ! - missing data increase rank in “transform domain” ! ! ‣ recovery using rank penalization scheme Saturday, November 11, 17

Nuclear-norm minimization convex relaxation of rank-minimization [Recht et. al., 2010] min || X || ∗ s.t. k A ( X ) � b k 2  ✏ X { sum of singular values of X [Rennie and Srebro 2005, Lee et. al. 2010, Recht and Re 2011] X = LR H n ry × n sy n f n ry × n sy n k n k R H n rx × n sx n = n f r x × n L X n f s x 18 X ∈ C n f × n rx × n sx × n ry × n sy L ∈ C n f × n rx × n sx × n k R ∈ C n f × n ry × n sy × n k Saturday, November 11, 17

[Rennie and Srebro 2005] Factorized formulation ! ‣ Upper-bound on nuclear norm is defined as ! 2 ! � ��  k LR H k ∗  1 L � � ! � � R 2 � � F ! where is sum of squares of all entries ! k . k 2 F ! ! ‣ choose explicitly & avoid costly SVD’s ! k ! Saturday, November 11, 17

Survey information –– coil acquisition 20 Saturday, November 11, 17

Old vs new Conventional acquisition random coil acquisition from https://www.slb.com/~/media/Files/resources/oilfield_review/ors08/aut08/shooting_seismic_surveys_in_circles.pdf 21 Saturday, November 11, 17

Acquisition mask – non-canonical matrix (10 x10 km) 22 Saturday, November 11, 17

Acquisition information 3D overthrust model, 5km x 12km x 12km 10404 sources @ 100m 40804 receivers @ 50m Time length : 3 seconds @ 0.004s Interpolation from 1-50 Hz Saturday, November 11, 17

Acquisition information 3D overthrust model, 5km x 12km x 12km 10404 sources @ 100m Unknown 20k X 20k matrix for each frequency! 40804 receivers @ 50m Time length : 3 seconds @ 0.004s Interpolation from 1-50 Hz Saturday, November 11, 17

Frequency slice @ 7Hz ground truth # 10 4 common source gather 0.2 0 0.4 2 0.6 4 Ry (km) 0.8 6 Sy-Ry 1 8 1.2 10 1.4 0 2 4 6 8 10 Rx (km) 1.6 1.8 2 0.5 1 1.5 2 24 # 10 4 Sx-Rx Saturday, November 11, 17

Frequency slice @ 7Hz observed # 10 4 common source gather 0.2 0 0.4 2 0.6 4 Ry (km) 0.8 6 Sy-Ry 1 8 1.2 10 1.4 0 2 4 6 8 10 Rx (km) 1.6 1.8 2 0.5 1 1.5 2 25 # 10 4 Sx-Rx Saturday, November 11, 17

Frequency slice @ 7Hz interpolated × 10 4 common source gather 0.2 0 0.4 2 0.6 4 Ry (km) 0.8 6 Sy-Ry 1 8 1.2 10 1.4 0 2 4 6 8 10 Rx (km) 1.6 1.8 2 0.5 1 1.5 2 26 × 10 4 Sx-Rx Saturday, November 11, 17

Frequency slice @ 7Hz residual × 10 4 common source gather 0.2 0 0.4 2 0.6 4 Ry (km) 0.8 6 Sy-Ry 1 8 1.2 10 1.4 0 2 4 6 8 10 Rx (km) 1.6 1.8 2 0.5 1 1.5 2 27 × 10 4 Sx-Rx Saturday, November 11, 17

Common source gather ground truth 0 2 Receiver-Y (km) 4 6 8 10 0 2 4 6 8 10 Receiver-X (km) 0 0 0.5 0.5 1 1 Time (s) Time (s) 1.5 1.5 2 2 2.5 2.5 3 3 28 0 2 4 6 8 10 0 2 4 6 8 10 Receiver-X (km) Receiver-Y (km) Saturday, November 11, 17