4 ProbabilisticSeismicHazardAnalysis(PSHA)reliesontwowidelyutilizedapproacheswithhighcomputationaldemands: 5 (a)Riemannsumand(b)conventionalMonteCarlo(MC)integration.Thefirstrequiressufficientlyfineslicesacrossmag- 6 nitude, distance, and ground motion, and the second requires extensive synthetic earthquake catalogs to compute seismic 7 hazardsaccurately.Theseapproachesarenotablyresource-intensiveforlow-probabilityseismichazards,e.g.,upto 108MC 8 samplesforahazardwith 10−4probabilitytoachievecoefficientofvariation(COV)of1%.Here,wepresentanovelframework 9 to compute hazard and deaggregation with unprecedented computational efficiency. We formulate Adaptive Importance 10 Sampling(AIS)PSHAtoapproximateoptimalimportantsampling(IS)distributionsanddramaticallyreducethesizeofsyn- 11 theticearthquakecatalogs(i.e.,numberofMCsamples)toestimatehazards.Weevaluatetheeffectivenessandreliabilityof 12 ourproposedmethodusingcomprehensivetestproblemsfromthePacificEarthquakeEngineeringResearchCenter(PEER) 13 forPSHAbenchmarks,encompassingvariousseismicsourcetypes,includingareal,fault,andcombinedones.Ourfindings 14 indicatethatthisnovelapproachsignificantlyoutpacesRiemannsumandtraditionalMCmethodswithcomputationsupto 15 >105and7.8×103timesfaster,respectively,whilemaintaininganstandarddeviationoftheestimatebelow2%.Moreover,we 16 showtheoreticallythatoptimalISdistributionsareequivalenttohazarddeaggregationdistributions.Empirically,weshow 17 ourapproximatedoptimalISandthedeaggregationdistributionsarecloselyalike,e.g.,withaKolmogorov–Smirnovstatistic 18 between0.017and0.113.WedevelopedourmethodologytohavebroadapplicabilityinPSHApractices,especiallyincases 19 requiringextensivecomputationalresourcestonavigatenumerouslogictreescenariosaddressingepistemicuncertainty. 20 1. Civil and Environmental Engineering, University of California, Berkeley, California, U.S.A, https://orcid.org/0000-0001-6873-5327 (SEH) https://orcid.org/0000- 0003-0322-7510 (LC) *Corresponding author: shoung@berkeley.edu Cite this article as Houng, Soung Eil and Luis, Ceferino (2022). Fast Probabilistic Seismic Hazard Analysis through Adaptive Importance Sampling, Bull. Seismol. Soc. Am.0, 1–35, doi: 00.0000/0000000000. © Seismological Society of America Volume 0 Number 0 XXXX XXXX www.bssaonline.org •1 KEY POINTS •We introduce an adaptive importance sampling algorithm for Probabilistic Seismic Hazard Analysis (PSHA). •The proposed algorithm is up to 103and 105times more efficient than conventional Monte Carlo and Riemann sum. •The algorithm also facilitates the straightforward implementation of hazard deaggregation. Supplemental Material INTRODUCTION 21 ProbabilisticSeismicHazardAnalysis(PSHA)hasbecomeafoundationalmethodfordeterminingseismicdesignlevelsand 22 conductingregionalseismicriskanalysessinceitsfirstinception(Cornell,1968;U.S.NuclearRegulatoryCommission,2007; 23 McGuire,2008;ASCE,2022).ThedevelopmentofPSHAwasdrivenbythenecessityforaprobabilisticframeworktoquantify 24 seismichazards,acknowledgingtheunpredictablenatureofearthquakelocation,magnitudes,andgroundmotionintensities 25 (Esteva,1967;Cornell,1968).InPSHA,earthquakelocation,magnitude,andgroundmotionsaretreatedasrandomvariables, 26 facilitating the computation of annual exceedance probabilities at various ground motion intensities. Integrating hazard 27 with fragility curves further enables to determine the annual probabilities of structural failures, thereby underscoring the 28 methodology’scriticalroleinriskassessment(Kennedyetal.,1980;Ceferinoetal.,2020;Silvaetal.,2020;Bakeretal.,2021; 29 PapadopoulosandBazzurro,2021;AroraandCeferino,2023). 30 Since we cannot solve PSHA analytically due to the complexity in seismic source and ground motion models, numerous 31 researchers have developed computer softwares for PSHA computation (Cornell, 1968; McGuire, 1976; Field et al., 2003; 32 Ordaz et al., 2013; Pagani et al., 2014). According to the PSHA computer code verification project (Thomas et al., 2010; 33 Haleetal.,2018),existingsoftwareprimarilyemploysRiemannsummationfornumericalintegrationofPSHA,whichpar- 34 titions the earthquake magnitude, location, and ground motion random variable into fine grids to approximate the actual 35 integration. The Riemann summation offers robust PSHA integration with sufficiently dense grids. However, this method 36 generallyincursasignificantcomputationalloadexponentiallyincreasingwiththenumberofgridsanddimensionsinmulti- 37 dimensionalintegrations(PhilippeandRobert,2001).Furthermore,theresultsarehighlysensitivetothechosengriddesign 38 (e.g.,theinitialpointofthegrid,gridspacing),leadingtosignificantdeviationsfromonesoftwaretoanother,especiallyfor 39 lowexceedanceprobabilities(Thomasetal.,2010;Haleetal.,2018). 40 Alternatively,othersoftwaresadoptedMonte-Carlo(MC)integrationforPSHA(AssatouriansandAtkinson,2013,2019). 41 MCintegrationcalculatesexceedanceprobabilitiesbygeneratingsyntheticearthquakecatalogsbasedontheseismicsource 42 andgroundmotionmodelstoevaluatetherecurrenceofvariousgroundmotionintensities(Musson,2000).MCintegration’s 43 primary advantage lies in its straightforward concept compared to Riemann summation, without the need to divide the 44 integration range into small slices (Musson, 2000; Dick et al., 2013). Nevertheless, MC framework requires a substantially- 45 2• www.bssaonline.org Volume 0 Number 0 Month Y ear long synthetic catalog to accurately estimate hazard from rare events, making the computation costly, especially for large 46 groundmotionswithlowexceedanceprobabilities,e.g., 𝑝<10−4/yr(Kroeseetal.,2014). 47 ImportanceSampling(IS)canofferasolutiontothisrareeventsimulation(TokdarandKass,2010).ISwasinitiallyintro- 48 ducedinstatisticalphysics(HammersleyandMorton,1954)toimprovethecomputationalefficiencyofrareeventsimulation 49 that would otherwise require a large sample size with conventional MC. IS relies on identifying an appropriate probability 50 distribution (“IS distribution”) to explore low-probability spaces effectively. Researchers use the IS distribution to sample 51 rare events with a higher likelihood than conventional MC and then correct their frequency through weights, significantly 52 reducingthenumberofsamplestocomputelowprobabilities(Robertetal.,1999).However,findingsuchanappropriatedis- 53 tribution can be challenging because there is no optimal sampling density that is universally applicable; rather, the proper 54 selectionofsamplingdensitydependsontheproblembeingsolved.Thus,manynumericalexperimentsareoftenconducted 55 firstthroughtrial-and-errortoidentifyISdistributions,whichisstillcomputationallyexpensive. 56 Inregionalseismicriskanalysis,numerousstudieshavebeenconductedtosamplehazard-consistentearthquakeground 57 motions(CrowleyandBommer,2006;Kiremidjianetal.,2007;JayaramandBaker,2010;HanandDavidson,2012;Manzour 58 etal.,2016;Christouetal.,2018;Kavvadaetal.,2022).Kiremidjianetal.(2007)firstintroducedISdistributionsthatsamples 59 large-magnitudeearthquakeswithahighprobabilitytoreducethecomputationalburdenofseismichazardandriskanalyses. 60 Jayaram and Baker (2010) expanded the approach by defining IS distributions to sample high-intensity ground motions. 61 However, Jayaram and Baker (2010) highlighted the computational challenges to identify an effective IS distribution and 62 ended up using K-mean clustering to reduce the number of ground motion samples. Rahimi and Mahsuli (2019) applied 63 system reliability methods to calculate PSHA. They selected the IS sampling density as the normal distribution centered at 64 the"designpoint"derivedfromthefirst-andsecond-orderreliabilitymethod. 65 To find effective IS distributions, computational statisticians have developed a general framework to find them through 66 iterativealgorithmsdenominated“AdaptiveImportanceSampling(AIS)”(Bugalloetal.,2017).AISalgorithmsleveragethe 67 fact that the optimal density is proportional to the integrand of the problem to find the sampling density which minimizes 68 the variance of the estimate. Various AIS algorithms have been introduced, e.g., cross-entropy based AIS, Vegas, Divonne, 69 andMiser(Lepage,1978;FriedmanandWright,1981;Rubinstein,1997;PressandFarrar,1990;RubinsteinandKroese,2004; 70 Bugalloetal.,2017).TheapplicationofAISiswidelyadoptedtosolvetheintegrationwithhighdimensions,suchasinthe 71 fieldofstatisticalphysics,finance,reliabilityengineering,andsignalprocessing(AuandBeck,2001;KappenandRuiz,2016; 72 NietoandRuiz,2016;Bugalloetal.,2017).Althoughpreviousstudiessuggesttheuseofimportancesampling(IS)forPSHA 73 calculation(JayaramandBaker,2010;RahimiandMahsuli,2019),noresearchhasbeenpublishedregardingtheuseofAIS 74 for this purpose. AIS can provide a general methodology for identifying an appropriate IS distribution for seismic hazards, 75 eliminating computationally expensive experiments in the regular IS approach, whose efficiency typically depends on the 76 researcher’sexperienceinthefield. 77 Volume 0 Number 0 Month Y ear www.bssaonline.org •3 ManyAISalgorithmshavebeenpublished(Lepage,1978;FriedmanandWright,1981;Rubinstein,1997;PressandFarrar, 78 1990;RubinsteinandKroese,2004;Bugalloetal.,2017).However,westudytheVEGASalgorithmforPSHAtoleverageits 79 straightforwardmathematicalframeworkandfastconvergence.Thus,thisstudyintroducesanovelcomputationalmethod 80 forPSHAcurvecalculationusingAdaptiveImportanceSampling(AIS)VEGASalgorithm(Lepage,1978,2021).Thispaper 81 alsoshowsthatAISfacilitateshazarddeaggregation,therelativecontributionofeachrandomvariablestotheoverallhazard 82 (Bazzurro and Cornell, 1999), because the optimal IS density resembles the contributions of each variable—magnitude, 83 distance,andgroundmotion—tothehazard.Thus,ourproposedframeworkalsoenhancesthecomputationalefficiencyof 84 deaggregationestimates,whichotherwisewouldrequireadditionalmemoryresourcesandcomplexityinthecomputercode 85 with traditional calculation methods. We explore three key aspects of the method: 1) the enhancement of computational 86 efficiencythattheAISalgorithmoffersovertraditionalmethods;2)theaccuracyoftheestimatesprovidedbythisapproach; 87 and3)theprocessofobtaininghazarddeaggregationthroughAIS.WepresentthetheoreticalbackgroundofAISPSHAand 88 validatethemethodthroughnumericalexamples. 89 MATHEMATICAL FORMULATION 90 Probabilistic Seismic Hazard Analysis (PSHA) 91 Atasiteofinterest,theannualfrequencyofgroundmotionexceedancefromasinglesourcecanbecalculatedas: 92 𝜆(𝑋>𝑎)=𝜈∫𝜀max 𝜀min∫𝑟max 𝑟min∫𝑚max 𝑚min𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑑𝑚𝑑𝑟𝑑𝜀 (1) , where𝜆(𝑋>𝑎)is annual rate that ground motion, 𝑋, exceeds the target ground motion intensity, 𝑎, e.g., peak ground 93 acceleration. 𝜈istheannualrateofearthquakeoccurrencegreaterthan 𝑚minfromthesource, 𝑚istheearthquakemagnitude, 94 𝑚minand𝑚maxareminimumandmaximummagnitudesconsideredforthesource, 𝑟isthesource-to-sitedistance, 𝑟minand 95 𝑟maxare minimum and maximum source-site distances, 𝜀is a standard normal random variable for generating earthquake 96 ground motion, 𝜀minand𝜀maxare minimum and maximum 𝜀(generally, 𝜀max≥6and𝜀min≤−6; Bommer and Abrahamson 97 (2006)),𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)is joint probability density function (PDF) of 𝑚,𝑟, and𝜀,𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)is indicator function that 98 takes 1 when 𝑋>𝑎, otherwise, 0. The ground motion 𝑋given𝑀,𝑅, and ℰis generally calculated using ground motion 99 models (Bozorgnia et al., 2014; Goulet et al., 2021). The models usually assume the log-normal distribution for the ground 100 motion given explanatory variables such as 𝑀and𝑅. Naturally, these models provide the mean and standard deviation of 101 logarithmicgroundmotion.Thus,therandomgroundmotioniscalculatedas: 102 log𝑋=𝜇(𝑀,𝑅)+ℰ𝜎(𝑀,𝑅) (2) ,where𝜇and𝜎aremeanandstandarddeviationoflogarithmicearthquakegroundmotion.Bytakingexponentialonboth 103 sidesof Eq.(2),thegroundmotion 𝑋canbecalculatedas 104 4• www.bssaonline.org Volume 0 Number 0 Month Y ear 𝑋=𝑒𝜇(𝑀,𝑅)+ℰ𝜎(𝑀,𝑅) Ifweassumethatthegroundmotionrandomvariable 𝜀isindependentwithrespecttothe 𝑚and𝑟(McGuire,1995),Eq.(1) 105 canbemodifiedas 106 𝜆(𝑋>𝑎)=𝜈∫𝜀max 𝜀min∫𝑟max 𝑟min∫𝑚max 𝑚min𝐼(𝑥>𝑎|𝑚,𝑟,𝑒)𝑓ℰ(𝜀)𝑓𝑀,𝑅(𝑚,𝑟)𝑑𝑚𝑑𝑟𝑑𝜀 =𝜈∫𝜀max 𝜀min∫𝑅max 𝑅min∫𝑀max 𝑀min𝐼(𝑥>𝑎|𝑚,𝑟,𝑒)𝑓ℰ(𝜀)𝑓𝑅|𝑀(𝑟|𝑚)𝑓𝑀(𝑚)𝑑𝑚𝑑𝑟𝑑𝜀 (3) ,where𝑓𝑀(𝑚),𝑓𝑅|𝑀(𝑟|𝑚),𝑓ℰ(𝜀)arePDFof𝑚,𝑟given𝑚,and𝜀.Underpointsourceassumption,thedistance 𝑟andmagnitude 107 𝑚becomeindependentrandomvariables.Thus,theseismichazardisgivenby 108 𝜆(𝑋>𝑎)=𝜈∫𝜀max 𝜀min∫𝑟max 𝑟min∫𝑚max 𝑚min𝐼(𝑥>𝑎|𝑚,𝑟,𝜀)𝑓ℰ(𝜀)𝑓𝑅(𝑟)𝑓𝑀(𝑚)𝑑𝑚𝑑𝑟𝑑𝜀 (4) Thetotalseismichazardfrommultipleseismicsources(e.g.,differentfaults)isthesumofeach.Thus, 109 Λ(𝑋>𝑎)=𝑛∑ 𝑖=1𝜆𝑖(𝑋>𝑎) (5) ,whereΛ(𝑋>𝑎)isthetotalannualfrequencyofexceedanceofgroundmotion, 𝑎,𝑖isindexforseismicsources,and 𝑛isthe 110 totalnumberofseismicsources.UndertheassumptionofPoissonprocess,theannualprobabilityofexceedance, 𝜙(𝑋>𝑎), 111 canbeconvertedfromannualfrequencyofexceedance( Λ;Eq.(5))as 112 𝜙(𝑋>𝑎)=1−𝑒−Λ(𝑋>𝑎)(6) Hazard deaggregation 113 Wecanalsodeaggregatethetotalhazard(Eq.(6))tobetterunderstandtheearthquakesthatcontributemosttothehazard. 114 Deaggregationisalsousedtodeveloptheselectseismicrecords(e.g.,fromtheearthquakesthatcontributemosttothehazard) 115 and conduct non-linear time-history analyses for the design of many critical buildings (Bazzurro and Cornell, 1999; U. S. 116 NuclearRegulatoryCommission,2007).Mathematically,deaggregationofthehazardisthejointprobabilitydistributionof 117 the𝑚,𝑟,and𝜀conditionalondifferentlevelsofhazards 𝑎toquantifythecontributionsofeachcomponent.Thedeaggregation 118 ofPSHAcanbeformulatedusingBayes’theoremas 119 Volume 0 Number 0 Month Y ear www.bssaonline.org •5 𝑃(𝑚,𝑟,𝜀|𝑋>𝑎)=𝑃(𝑋>𝑎∩𝑚,𝑟,𝜀) 𝑃(𝑋>𝑎) =𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀) 𝑃(𝑋>𝑎) =𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀)∑𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀) =𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀)∑𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀)(7) ,where𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)istheprobabilityofgroundmotion 𝑋isgreaterthan 𝑎given𝑚,𝑟,and𝜀,𝑃(𝑚,𝑟,𝜀)isjointprobability 120 of𝑚,𝑟, and𝜀, and𝑃(𝑋>𝑎)is the total probability that the ground motion is greater than 𝑎, which is is the summation of 121 𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)𝑃(𝑚,𝑟,𝜀)over all𝑚,𝑟, and𝜀. Note that 𝑃(𝑋>𝑎|𝑚,𝑟,𝜀)can be expressed as an indicator function, 𝐼(𝑋> 122 𝑎|𝑚,𝑟,𝜀),becausetheprobabilityofgroundmotion 𝑋greaterthan 𝑎canonlybe1or0given 𝑚,𝑟,and𝜀(seeEq.(2)). 123 By replacing probability mass function, 𝑃(𝑚,𝑟,𝜀), with probability density function, 𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀), and change the 124 summationintointegration,Eq.(7)canbeexpressedas: 125 𝑓(𝑚,𝑟,𝜀|𝑋>𝑎)=𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) ∭𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑑𝑚𝑑𝑟𝑑𝜀(8) ByEq.(1),thedenominatorof Eq.(8)equals 𝜆∕𝜈.Thus, 126 𝑓(𝑚,𝑟,𝜀|𝑋>𝑎)=𝜈⋅𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) 𝜆(9) Thisequationshowsthatthecontributionofspecific 𝑚,𝑟,and𝜀canberepresentedbytheratioofthepartialsumofthe 127 given𝑚,𝑟, and𝜀to the total hazard. Note that 𝜆and𝜈are constant. Thus, the hazard deaggregation, 𝑓(𝑚,𝑟,𝜀|𝑥>𝑎), is 128 proportionalto 𝐼(𝑥>𝑎|𝑚,𝑟,𝜀)𝑓(𝑚,𝑟,𝜀): 129 𝑓(𝑚,𝑟,𝜀|𝑥>𝑎) ∝𝐼(𝑥>𝑎|𝑚,𝑟,𝜀)𝑓(𝑚,𝑟,𝜀) CURRENT NUMERICAL SOLUTIONS 130 The Riemann summation ThismethodcomputesPSHAcurvesbysummingtheareasofpartitioned( 𝑚,𝑟,𝜀)cuboids.The 131 RiemannsummationforEq.(1)canbeexpressedas: 132 𝜆(𝑋>𝑎)=𝜈𝑁𝜀∑ 𝑘=1𝑁𝑟∑ 𝑗=1𝑁𝑚∑ 𝑖=1𝐼(𝑥>𝑎|𝑚𝑖,𝑟𝑗,𝜀𝑘)𝑓𝑀,𝑅,ℰ(𝑚𝑖,𝑟𝑗,𝜀𝑘)∆𝑚∆𝑟∆𝜀 (10) 6• www.bssaonline.org Volume 0 Number 0 Month Y ear TABLE 1 Comparison of time complexity of various PSHA algorithms Algorithm 𝑁𝑚 𝑁𝑟 𝑁𝜀 𝑁𝑠 𝑁𝑎 Riemann sum 𝑁 𝑁 𝑁 - 1 Conventional MC 1 1 1 𝑁 1 IS MC 1 1 1 𝑁 1 VEGAS AIS (this study) 1 1 1 𝑁 𝑁 𝑁𝑚,thenumberofmagnitudegrids; 𝑁𝑟,thenumberofdistancegrids; 𝑁𝜀,thenumberof 𝜀grids;𝑁𝑠,thenumberof (𝑚,𝑟,𝜀) samples;𝑁𝑎,thenumberofgroundmotionintensityofinterest(e.g.,0.1g) , where∆𝑚,∆𝑟, and∆𝜀are grid step size for Riemann sum, 𝑁𝑚,𝑁𝑟, and𝑁𝜀are the total number of grids, satisfying 𝑁𝑥∆𝑥 133 =𝑥max-𝑥𝑚𝑖𝑛, where𝑥= {𝑚,𝑟,𝜀}. Note that 𝑓𝑀,𝑅,ℰ(𝑚𝑖,𝑟𝑗,𝜀𝑘)∆𝑚∆𝑟∆𝜀is equivalent to the probability at ( 𝑚,𝑟,𝜀) such that 134 |𝑚−𝑚𝑖|<∆𝑚∕2,|𝑟−𝑟𝑗|<∆𝑟∕2,and|𝜀−𝜀𝑘|<∆𝜀∕2.TheaccuracyoftheRiemannsummationdependsonthegridsize. 135 Utilizingfinergridsenhancestheaccuracyofthesummation.However,thecomputationtimeisinverselyproportionaltothe 136 gridstepsizeandthusproportionaltothenumberofgrids.Forthree-dimensionalPSHAsummation,thecomputationtime 137 scales with 𝑁𝑚×𝑁𝑟×𝑁𝜀(Table 1). Notably, since the distance PDF ( 𝑓𝑅(𝑟)) cannot be analytically determined in practice, 138 integration often extends over latitude ( 𝜙), longitude ( 𝜓), and depth ( 𝑧), increasing the dimensions from three ( 𝑚,𝑟,𝜀) to 139 five (𝑚,𝜙,𝜓,𝑧,𝜀). Therefore, Riemann summation for seismic hazard becomes even more computationally intensive, a 140 phenomenonknownas"thecurseofdimensionality"duetotheexponentialincreaseincomputationtimewiththenumber 141 ofdimensions(NovakandRitter,1997). 142 Volume 0 Number 0 Month Y ear www.bssaonline.org •7 Conventional Monte-Carlo (MC) ThismethodsimulatesmanysyntheticearthquakegroundmotionsandcalculatesPSHA 143 byassessingthefrequencywithwhichgroundmotionintensitiesexceedacertainthreshold.MCPSHAiscomputedas 144 ̂𝜆(𝑋>𝑎)=𝜈 𝑁𝑁∑ 𝑖=1𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖) =1 𝑇𝑁∑ 𝑖=1𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖) (11) ,where𝑋𝑖denotesthesimulatedgroundmotions,and 𝑁isthetotalnumberofsamples. 𝑀𝑖,𝑅𝑖,ℰ𝑖arerandomsamplesfrom 145 𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀),𝑇representstheequivalentcatalogdurationequalto 𝑁∕𝜈.MCPSHAisunbiasedbecausetheexpectationof 146 ̂𝜆(𝑋>𝑎)isthesameas 𝜆: 147 𝐸[̂𝜆]=𝐸⎡ ⎢ ⎣𝜈 𝑁𝑁∑ 𝑖=1𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖)⎤ ⎥ ⎦ =𝜈 𝑁𝑁∑ 𝑖=1𝐸[𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖)] =1 𝑁𝑁∑ 𝑖=1𝜈∫𝐼(𝑋𝑖>𝑎|𝑚𝑖,𝑟𝑖,𝜀𝑖)𝑓𝑀,𝑅,ℰ(𝑚𝑖,𝑟𝑖,𝜀𝑖)𝑑𝑚𝑖𝑑𝑟𝑖𝑑𝜀𝑖 =1 𝑁𝑁∑ 𝑖=1𝜆(𝑋>𝑎) =𝜆(𝑋>𝑎) (12) Thevarianceof ̂𝜆,thevariabilityofeachMCestimates,canbeobtainedasfollows(AppendixA): 148 VAR[̂𝜆]=𝜈𝜆−𝜆2 𝑁(13) NotethatVAR [̂𝜆]isalwayspositivesince 𝜈≥𝜆.Then,thestandarddeviationoftheestimate, ̂𝜆,canbeobtainedas: 149 𝜎[̂𝜆]=√ 𝜈𝜆−𝜆2 𝑁(14) , which provides the absolute uncertainty about MC PSHA estimates. However, we do not want to fully rely on 𝜎[̂𝜆]to 150 comparedifferentexceedanceprobabilities.Forexample,supposeweareinterestedintwodifferentexceedancefrequencies, 151 𝜆1=10−1/yr and𝜆2=10−2/yr when𝜈= 1/yr and 𝑁= 100. Then, the variance (uncertainty) of the two MC estimates, ̂𝜆1152 and̂𝜆2, are𝜎1=3×10−2and𝜎2∼1×10−2, respectively (Eq. (14)). Here, someone might argue that uncertainty of ̂𝜆1MC 153 estimateislargerthanthatof ̂𝜆2because𝜎1isgreaterthan 𝜎2.However,theuncertaintyof 𝜆1shouldbeconsideredsmaller 154 thanthatof 𝜆2consideringthetargettruevalueofeachestimate.Thatis,0.1 ±0.03(̂𝜆1)isbetterestimatethan0.01 ±0.01 155 (̂𝜆2).Instead,itisbettertousethecoefficientofvariation,COV,thestandarddeviationnormalizedbyitsmean 156 8• www.bssaonline.org Volume 0 Number 0 Month Y ear 𝐶𝑂𝑉=𝜎[̂𝜆] 𝐸[̂𝜆]=√ 𝜈−𝜆 𝑁𝜆(15) Smaller COV indicates that the PSHA estimate is more accurate relative to its actual value, e.g., 5 %COV means that the 157 hazardestimatesarewithin ±5%ofitstruevaluewitha68%probability,assumingnormallydistributedMCestimatesasper 158 thecentrallimittheorem.COVisinverselyproportionaltothesquarerootofthenumberofMCsamples( 𝑁),indicatingthat 159 increasing𝑁naturallyimprovestheestimate’saccuracy. 160 Thetargetexceedancefrequency, 𝜆,alsoaffectstheMCaccuracy.Thelower 𝜆leadstopooraccuracywithfixed 𝑁and𝜈.It 161 isintuitivelyreasonablethatsufficientlylongearthquakecatalogisrequiredtoaccuratelyestimatetheeventfromlongreturn 162 period. It is also notable that the variance of MC estimate is proportional to the square root of earthquake occurrence rate, 163 𝜈, indicating that MC PSHA is more challenging task in region with higher earthquake activity. This can also be explained 164 intuitively.Regionswithhighseismicactivityexperiencemoreearthquakesperyearthanlessactiveareas,therefore,when 165 calculatingthegroundmotionexceedanceprobabilityoverthesameperiod,moreearthquakesneedtobeconsideredinthese 166 activeregions. 167 AnadvantageofMCPSHAisthatthecomputationaltimeisdependentonthenumberofsamples 𝑁,circumventingthe 168 inherentdimensionalityprobleminRiemannsummation(Table1).Inotherwords,thecomputationtimeofMCsimulation 169 is independent of the number of grids ( 𝑁𝑚,𝑁𝑟, and𝑁𝜀in Eq. (10)). Therefore, we could utilize fine joint probability mass 170 functionsformoreprecisehazardestimationwithoutanincreaseincomputationalburden.Ifclosed-formprobabilitydensity 171 functionsareavailable,PSHAcanbeimplementedwithoutapproximateddiscretization.Also,thecomputationtimeofMCis 172 independentofthenumberofgroundmotionsofinterest, 𝑁𝑎(Table1).Wecanutilizethegeneratedsyntheticgroundmotion 173 catalogtotallytheexceedanceeventsforallthegroundmotionsofinterest,thoughwestillshouldrepeatthecomputation 174 alongwiththenumberofMCsamples. 175 ConventionalMCfacesextremecomputationalchallengesforlowprobabilitiesbecausethenumberofsamplesrequired 176 toachievelowCOVsdramaticallyincreases.WecancomputetherequirednumberofsamplesbyrearrangingEq.(15): 177 𝑁=1 (𝐶𝑂𝑉)2×𝜈−𝜆 𝜆(16) Forsmall𝜆,i.e.,𝜈>>𝜆,theEq.(16)canbeapproximatedas: 178 𝑁∼1 (𝐶𝑂𝑉)2×𝜈 𝜆(17) Volume 0 Number 0 Month Y ear www.bssaonline.org •9 To demonstrate conventional MC’s extreme computational demands, we can calculate 𝑁for a target annual exceedance 179 frequency,𝜆,of10−4peryearandaannualearthquakerate 𝜈=1/yr,whichistypicalvaluesinPSHApractice(Coppersmith 180 etal.,2014).FromEq.(17),weneed 𝑁∼108toachieveCOV=1%. 181 Importance Sampling (IS) Integration 182 IS is a generalization of MC theory. Consider a random variable 𝑥that follows a probability function, 𝑓𝑋(𝑥). The expected 183 valueofafunction 𝑢(𝑥),denotedas 𝑆,isdefinedby 184 𝑆=∫𝑢(𝑥)𝑓𝑋(𝑥)𝑑𝑥 (18) TheMCestimateis 185 ̂𝑆=1 𝑁𝑁∑ 𝑖=1𝑢(𝑋𝑖) (19) ,where𝑋𝑖issampledfrom 𝑓𝑋(𝑥). 186 Byintroducinganarbitraryprobabilityfunction, 𝑞𝑋(𝑥),Eq.(18)canbeequivalentlyexpressedas 187 𝑆=∫𝑢(𝑥)𝑓𝑋(𝑥) 𝑞𝑋(𝑥)𝑞𝑋(𝑥)𝑑𝑥 (20) WerestricttheintegrationrangeinEq.(20)where 𝑓𝑋(𝑥)≠0because𝑥suchthat𝑓𝑋(𝑥)=0doesnotcontributetotheinte- 188 gration. Then, 𝑞𝑋(𝑥)can be any distribution with nonzero density in the integration range. Eq. (20) provides important 189 implications in MC estimation. We can get the solution to Eq. (18) by estimating the expected value of 𝑢(𝑥)𝑓𝑋(𝑥)∕𝑞𝑋(𝑥) 190 where𝑥follows the distribution 𝑞𝑋(𝑥). This approach is highly useful in numerical integration, especially when sampling 191 from𝑓𝑋(𝑥)is challenging or the population of 𝑓𝑋(𝑥)is extremely low in the region of importance, e.g., the exceedance of 192 PGAgreaterthan1gismostlycontributedbythegroundmotionsamplesfromlargemagnitude( 𝑚)earthquakesoccurred 193 atclosedistance( 𝑟)withlarge 𝜀,whichalltypicallycorrespondtothelowprobabilityregion. 194 TheISMCestimateis: 195 ̂𝑆=1 𝑁𝑁∑ 𝑖=1𝑢(𝑋𝑖)𝑓𝑋(𝑋𝑖) 𝑞𝑋(𝑋𝑖)(21) , where𝑋𝑖s are sampled from 𝑞𝑋(𝑥), the proposed (or new) IS sampling density function. If 𝑞𝑋(𝑥)is equal to the original 196 distribution, 𝑓𝑋(𝑥), Eq. (21) simplifies to the conventional Monte-Carlo (Eq. (19)). The ratio 𝑓𝑋(𝑥)∕𝑞𝑋(𝑥), known as the 197 importanceweight( 𝑤𝑖),adjustsforthechangeinsamplingdistribution. 198 IS PSHA Researchers have applied IS to PSHA using different sampling functions (Jayaram and Baker, 2010; Rahimi and 199 Mahsuli,2019).ForISPSHA,Eq.(1)canbereformulatedbyintroducinganewsamplingjointdensityfunction, 𝑞𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) 200 10 • www.bssaonline.org Volume 0 Number 0 Month Y ear asfollows 201 𝜆(𝑋>𝑎)=𝜈∫𝜀max 𝜀min∫𝑟max 𝑟min∫𝑚max 𝑚min𝐼(𝑥>𝑎|𝑚,𝑟,𝑒)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) 𝑞𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑞𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑑𝑚𝑑𝑟𝑑𝜀 (22) TheISMCestimatorofequation(22)is 202 ̂𝜆(𝑋>𝑎)=𝜈 𝑁𝑁∑ 𝑖=1𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,𝐸𝑖)𝑓𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,𝐸𝑖) 𝑞𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,𝐸𝑖) ThemeanofISestimator ̂𝜆isalsounbiasedlikeconventionalMCbecause 203 𝐸𝑞[̂𝜆]=𝜈 𝑁𝑁∑ 𝑖=1𝐸𝑞[𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖)𝑓𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖) 𝑞𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖)] =1 𝑁𝑁∑ 𝑖=1𝜈∫𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) 𝑞𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑞𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)𝑑𝑚𝑑𝑟𝑑𝜀 =1 𝑁𝑁∑ 𝑖=1𝜆 =𝜆 (23) Also,thevarianceoftheISPSHAestimateswithrespecttothetrue 𝜆isgivenby 204 VAR[̂𝜆]=1 𝑁⎛ ⎜ ⎝𝜈2𝐸⎡ ⎢ ⎣(𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖)𝑓𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖) 𝑞𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖))2⎤ ⎥ ⎦−𝜆2⎞ ⎟ ⎠(24) Hence,theCOVofISPSHAestimatecanbeexpressedas 205 𝐶𝑂𝑉=√ √√√√√√𝜈2𝐸[(𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,ℰ𝑖)𝑓𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖) 𝑞𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,ℰ𝑖))2 ]−𝜆2 𝑁𝜆2(25) Equivalency of optimal IS density and hazard deaggregation 206 FromEq.(24),wecanspecificallychooseanewsamplingdensity 𝑞∗ 𝑀,𝑅,ℰthatmakesVAR[ ̂𝜆]=0, 207 𝑞∗ 𝑀,𝑅,ℰ=𝜈𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,𝐸𝑖)𝑓𝑀,𝑅,ℰ(𝑀𝑖,𝑅𝑖,𝐸𝑖) 𝜆(26) Using𝑞∗ 𝑀,𝑅,ℰ, we could compute the true hazard, 𝜆, with only one MC sample because IS MC is unbiased (Eq. (23)) and 208 thevarianceiszero.Wecall 𝑞∗theoptimalISdensityforPSHAcalculation. 209 Volume 0 Number 0 Month Y ear www.bssaonline.org •11 We note that Eq. (26) is exactly the same as Eq. (9), indicating that the optimal density, 𝑞∗, is identical to the hazard 210 deaggregation. That is, if we find 𝑞∗, we are able to not only dramatically enhance the computational efficiency of seismic 211 hazardestimationbutalsoobtainthehazarddeaggregationdistributionsasaby-product. 212 Infact,wecanalsoseethisprofoundrelationshipbetweenthehazardanddeaggregationestimatesbyrearrangingEq.(9): 213 𝜆=𝜈⋅𝐼(𝑋>𝑎|𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) 𝑓(𝑚,𝑟,𝜀|𝑋>𝑎) Becausethisidentityholdsforanyvaluesof (𝑚,𝑟,𝜀),wecanobtainthehazardatgroundmotionlevel, 𝑎,withany(𝑚,𝑟,𝜀) 214 tripletifweknow 𝑓(𝑚,𝑟,𝜀|𝑋>𝑎).Theterm𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀)∕𝑓(𝑚,𝑟,𝜀|𝑋>𝑎)canbeinterpretedastheimportanceweightof 215 IS,and𝑓(𝑚,𝑟,𝜀|𝑋>𝑎),thehazarddeaggregation,istheoptimaldensity. 216 In existing numerical methods, the deaggregation can be more computationally expensive than the hazard. In Riemann 217 sum,wemustsaveeachsumelementinmemoryandallocatethoseelementsintoappropriatedeaggregationbins.Also,in 218 conventionalMC,weshouldsavethelongsyntheticgroundmotioncatalogwiththecorresponding( 𝑚,𝑟,𝜀)triplettoallocate 219 thoseintotheproperbins.Theseoperationsnecessitatesignificantcomputationalmemoryandtime.Therefore,theproperty 220 that the optimal density resembles the hazard deaggregation can be considered a huge benefit for hazard analysts. IS MC 221 with different densities other than 𝑞∗could also improve the computational efficiency compared to the conventional MC; 222 however,itcannotgiveanyinformationonthehazarddeaggregation(RahimiandMahsuli,2019;JayaramandBaker,2010). 223 Thoughtheuseofoptimaldensityisamajorbenefitincomputationofseismichazardandhazarddeaggregation,obtaining 224 itisnottrivialbecause 𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,𝐸𝑖)and𝜆inEq.(26)areunknownsbeforecalculationofPSHA.Thus,weproposea 225 newPSHAcomputationmethodtofindit. 226 ADAPTIVE IMPORTANCE SAMPLING PSHA 227 In“adaptive”importancesampling(AIS),weiterativelytraintheISdensitytofindtheoptimalone( 𝑞∗)byexploringimpor- 228 tantregionstocompute 𝐼(𝑋𝑖>𝑎|𝑀𝑖,𝑅𝑖,𝐸𝑖)and𝜆withareducednumberofMCsamples.AISmustbalancedifferentfactors 229 indetermininghowmanysamplesshouldbeusedtotrainoptimalISdensity.FewerMCsamplescanmakeAISfailasthey 230 will not allow us to explore the important regions effectively. On the other hand, many samples can impose large compu- 231 tational demands, even larger than those from conventional MC. In AIS, we must also consider computational costs are 232 proportionaltothenumberofiterationsforconvergence.Thus,itisimportanttoselectappropriatealgorithmsthatconverge 233 fast. 234 VEGAS Formulation for PSHA VEGASisanon-parametricadaptiveimportancesampling(AIS)algorithmthatiteratively 235 identifiestheoptimalproposaldensity, 𝑞∗(Lepage,1978,2021).Thealgorithmhasbeendevelopedandwidelyusedincom- 236 12 • www.bssaonline.org Volume 0 Number 0 Month Y ear putationalphysics(KersevanandRichter-Was,2013;Alwalletal.,2014),andiscurrentlyappliedtochemistry,astrophysics, 237 finance, and medical statistics (Campolieti and Makarov, 2007; Garberoglio and Harvey, 2011; Ray et al., 2011; Sanders, 238 2014). 239 The VEGAS algorithm is conceptually straightforward and is recognized for its rapid convergence, especially when the 240 randomvariablesinvolvedareindependent(Lepage,1978,2021).InPSHA,thevariable 𝜀isalwaysconsideredanindependent 241 variable(Eq.(3)).Additionally,thevariables 𝑚and𝑟arealsotreatedasindependentunderpointsourceassumption(Eq.(4)). 242 Notice, however, that when the finite-fault rupture model, in which the rupture dimension changes with magnitude, is 243 adopted,𝑚and𝑟canbecorrelated,andthedistributionof 𝑟isconditionalonthemagnitude 𝑚.Also,notethattheestimate 244 remainsunbiasedinthiscase,becauseitscalculationisstillwithintheISframework(Eq.(23)). 245 In three-dimensional integration, which is the case of PSHA, VEGAS employs 𝑁3cuboids that are independently parti- 246 tioned. The probability assigned to each cuboid and the total number of partitioned cuboids ( 𝑁3) are preserved across the 247 iteration steps. However, the IS sampling density changes because the algorithm updates the cuboid’s size depending on 248 its contribution to the integration. If a cuboid’s contribution is low, its size grows in the next step, lowering its probability 249 density.Conversely,whenacuboid’scontributionishigh,itshrinksinthefollowingstep,elevatingtheprobabilitydensity. 250 Ideally,wheneverycuboid’scontributiontotheintegrationbecomesidentical,wefindtheproposedoptimaldensity, 𝑞∗,and 251 thealgorithmisterminated.TheframeworktofindoptimalISdensitiesforPSHAusingVEGASalgorithmisexplainedwith 252 asimplepointseismicsourceexampleinthefollowingparagraphs. 253 First, we adopt an independently distributed joint probability function as IS density. Note that we could also include 254 correlations in the IS density, but such approach would increase computational memory and time demands, e.g., the com- 255 putationalcomplexityincreasesexponentiallywitheachaddeddimension,i.e., 𝑂(𝑁𝑑)(Lepage,1978).Incontrast,whenwe 256 assume independence, the computational complexity grows linearly with the number of dimensions, i.e., 𝑂(𝑁𝑑), making 257 multi-dimensionalintegrationinPSHAexceptionallyefficient.Thus, 258 𝑞(𝑚,𝑟,𝜀)=𝑞𝑀(𝑚)𝑞𝑅(𝑟)𝑞ℰ(𝜀) (27) Then,theintegrationrangesforeachvariable— 𝑚,𝑟,and𝜀—aredividedinto 𝑁gridswiththesamevolume.Thisdivision 259 isdesignedtogeneratecuboidsofconstantprobability: 260 𝑀∶𝑚𝑖−1≤𝑚<𝑚𝑖(𝑖=0,1,2,..,𝑁),∆𝑚𝑖=𝑚𝑖−𝑚𝑖−1 𝑅∶𝑟𝑗−1≤𝑟<𝑟𝑗(𝑗=0,1,2,..,𝑁),∆𝑟𝑗=𝑟𝑗−𝑟𝑗−1 ℰ∶𝜀𝑘−1≤𝜀<𝜀𝑘(𝑘=0,1,2,..,𝑁),∆𝜀𝑘=𝜀𝑘−𝜀𝑘−1(28) Volume 0 Number 0 Month Y ear www.bssaonline.org •13 Figure 1. VEGAS iterations of IS density in our AIS PSHA framework. Example for PSHA at PGA = 0.5 g when the site is located 10 km away from a point seismic source ( 𝑚min=5.0,𝑚max=8.0,𝑏-value = 1.0, 𝜈= 1.0/yr, GMM = Sadigh et al. (1997)). The total number of MC samples is 2,000, and the number of grids per axis is 10, constituting total 100 rectangles (cuboids in actual three dimensional PSHA integration). The vertical and horizontal black solid lines are boundaries of the 𝑚and𝜀rectangles, and gray dots are the MC samples. Starting from the initial same-size 𝑚and𝜀rectangles, their sizes are adjusted depending on the the contribution of each rectangle to the hazard. The final proposed structure of the rectangles gives a highly concentrated probability density at𝑚∼5and𝜀∼2. Note that the hazard estimates using initialized density, density after the first iteration, and the final proposed density are 0.0417, 0.0413, and 0.0378, while the true solution is 0.0385 The number of grids, 𝑁, is chosen to be 50 as suggested by Lepage (1978) and because we consider this value makes grids 261 sufficientlyfinetocapturetheactualdistributionoftheoptimaldensity 𝑞∗.TheprobabilityforeachcuboidfromEq.(28)is 262 set to1∕𝑁3to ensure that the initial IS density function is uniformly distributed across the entire domain. As a result, the 263 initialprobabilitydensity 𝑞(0)foraspecificcuboidis 264 𝑞(0)(𝑚𝑖−1≤𝑚<𝑚𝑖, 𝑟𝑗−1≤𝑟<𝑟𝑗, 𝜀𝑘−1≤𝜀<𝜀𝑘)=𝑞(0) 𝑀(𝑚)𝑞(0) 𝑅(𝑟)𝑞(0) ℰ(𝜀) =1 𝑁∆𝑚𝑖1 𝑁∆𝑟𝑗1 𝑁∆𝜀𝑘 Figure 1 illustrates an example of initial partitioning when 𝑁=10. The gray dots in the figure shows an example of MC 265 samplesfromthisinitial 𝑞0.Notethatforeasiervisualizationandunderstanding,Figure1illustratesanexamplein 𝑚-𝜀,the 266 two-dimensionalspace,andnotthethree-dimensional 𝑚,𝑟,and𝜀space. 267 Weupdate𝑞usingMCsamplestoultimatelymakeitconvergeto 𝑞∗(Eq.(26)).Asmentionedearlier,thesizeofeachcuboid 268 issubjecttochangewhiletheprobabilityofeachcuboidremainsconstant,i.e.,probabilitydensitychanges.Thisadjustment 269 is facilitated through a "subdivision-and-restoration" process (Lepage, 1978). In this process, 𝑖th grid is subdivided into 𝑛𝑖270 sub-grids, with 𝑛𝑖being proportional to the 𝑖th grid’s contribution to the overall integration, and restored to the original 271 number,𝑁, by merging 𝑁subgrid∕𝑁consecutive subgrids, where 𝑁subgridis the total number of sub-grids. Thus, the number 272 14 • www.bssaonline.org Volume 0 Number 0 Month Y ear ofsubdivisionsat 𝑖thgridis 273 𝑛𝑖=𝑁subgrid×⎛ ⎜ ⎝𝐻𝑖∆𝑚𝑖∑ 𝑖𝐻𝑖∆𝑚𝑖⎞ ⎟ ⎠,(𝑖=1,2,...,𝑁) (29) 274 𝑛𝑗=𝑁subgrid×⎛ ⎜ ⎝𝐻𝑗∆𝑟𝑗 ∑ 𝑗𝐻𝑗∆𝑟𝑗⎞ ⎟ ⎠,(𝑗=1,2,...,𝑁) (30) 275 𝑛𝑘=𝑁subgrid×⎛ ⎜ ⎝𝐻𝑘∆𝜀𝑘∑ 𝑘𝐻𝑘∆𝜀𝑘⎞ ⎟ ⎠,(𝑘=1,2,...,𝑁) (31) 𝑁subgridshouldbesufficientlylargerthan 𝑁toiteratetheISdensityeffectively,especiallywhenthegrids’contributionsto 276 thehazardfromthepreviousstagearehighlyheterogeneous.Inthisstudy,wechose 𝑁subgridtobe10,000,whichis200times 277 greaterthan 𝑁(=50).Thesecondtermof Eq.(29)-Eq.(31)representstheportionofeachgrid’scontributiontothehazard. 278 AccordingtoLepage(1978), 𝐻𝑖inEq.(29)is 279 𝐻𝑖=√ √√√∑ 𝑚𝑖−1<𝑀<𝑚𝑖𝐻2(𝑀,𝑅, ℰ;𝑎) 𝑞𝑟(𝑅)𝑞𝜀(ℰ)(32) where 280 𝐻(𝑀,𝑅, ℰ;𝑎)=𝐼(𝑋>𝑎|𝑀,𝑅, ℰ)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀) Note that𝐻(𝑀,𝑅, ℰ;𝑎)is the integrand of PSHA. Intuitively, 𝐻𝑖can be considered as the marginalization of the overall 281 contributionofthemagnitudedimensionwithinthe 𝑖thgrid.Notethatthesummationof 𝐻2sinEq.(32)isdoneoverallthe 282 sampleswhichfallwithin 𝑚𝑖−1and𝑚𝑖andthedivisionby 𝑞𝑟𝑞𝜀canbeconsideredanadjustmentforunevenlydistributed 𝑟 283 and𝜀samplestopurelycapturethemagnitudecontribution. 284 Similarly,𝐻𝑗and𝐻𝑘inEq.(30)andEq.(31)canbecalculatedas 285 𝐻𝑗=√ √√√∑ 𝑟𝑖−1<𝑅<𝑟𝑖𝐻2(𝑀,𝑅, ℰ;𝑎) 𝑞𝑚(𝑀)𝑞𝜀(ℰ) 286 𝐻𝑘=√ √√√∑ 𝜀𝑖−1<ℰ<𝜀𝑖𝐻2(𝑀,𝑅, ℰ;𝑎) 𝑞𝑟(𝑅)𝑞𝑚(𝑀) Intheearlyiterations,researchershavenoticednumericalinstabilitiesduetotherelativelypoorinformationontheinte- 287 grand(Lepage,1978,2021).However,researchershavealsofoundeffectivewaystomitigateitthroughsmoothing(Eq.(33)) 288 Volume 0 Number 0 Month Y ear www.bssaonline.org •15 anddamping(Eq.(34)). 289 𝑑𝑖∶=1∑𝑑𝑖×⎧ ⎪⎪ ⎨⎪⎪ ⎩(7𝑑𝑖+𝑑𝑖+1)∕8 if𝑖=0 (𝑑𝑖−1+6𝑑𝑖+𝑑𝑖+1)∕8if𝑖=1,2,...,𝑁−1 (𝑑𝑖−1+7𝑑𝑖)∕8 if𝑖=𝑁(33) ,where𝑑𝑖is𝐻𝑖∆𝑥𝑖∕∑ 𝑖𝐻𝑖∆𝑥𝑖,and𝑥are𝑚,𝑟,or𝜀inEq.(29)-(31). 290 𝑑𝑖∶=(1−𝑑𝑖 ln(1∕𝑑𝑖))𝛼 (34) ,where𝛼islearningrate.Weused 𝛼as1.0assuggestedbyLepage(2021).Throughnumericalexperimenting,weobserved 291 thisvalueallowedmostofthePSHAintegrationtoconvergewithinthreeiterations(asdescribedlater)withoutcausingany 292 numericalinstabilities. 293 Finally, based on calculated 𝑛𝑖,𝑛𝑗and𝑛𝑘(Eq. (29)-Eq. (31)), the number of grid is restored to the original size, 𝑁, by 294 merging𝑁subgrid∕𝑁consecutive subgrids. The restored grid is the updated density, 𝑞(1). Figure 1 illustrates this iteration. 295 The sizes of the grids are shrunk at 5<𝑚<6and𝜀∼2, indicating that the high contributions to the hazard on this range 296 (Eq.(29)-Eq.(31))incontrasttootherlessimportantregions,e.g.,at 𝑚>6.5and𝜀<1.Notethattheprobabilityofeachgrid 297 ispreservedsothatthenumberofMCsamplesinsideeachgridisalmostthesameregardlessofthegridsize. 298 Theiterativeprocessexplainedaboveisrepeateduntilnofurtherimprovementisobservedinthevarianceofthehazard 299 estimator(Eq.(24)).AnexampleofthegridstructureofthefinalproposedISdensityisshowninFigure1.Thefinalproposed 300 densityexhibitsfeatureswhichisexpectedbasedonintuition.First,thelowcontributionof 𝜀lessthan0isunderstandable. 301 Thisisbecausethesimulatedgroundmotionintensity,usingSadighetal.(1997),neverexceedsthetargetgroundmotionof 302 0.5𝑔atanymagnitude,evenat 𝑚=8.0,when𝜀islessthan0.1.Inaddition,thestrongcontributionof 𝑚∼5and𝜀∼2.2to 303 theintegrationmakessensebecausethesimulatedgroundmotionintensityexceedsthetargetgroundmotionof0.5gwhen 304 𝑚and𝜀reach5and2.2,respectively.Thedecreasingtrendofcontributionbeyond 𝑚∼5and𝜀∼2.2canbeinterpretedasthe 305 exponentialdecayintheprobabilityof 𝑚and𝜀.ThisresultforasimplePSHAexampleimpliesthattheproposedalgorithm 306 canbeutilizedformorecomplexPSHAintegration. 307 NUMERICAL EXAMPLES 308 WetestedourproposedAISPSHAonthecomprehensivebenchmarkproblemsets1.11and2.1fromthePEERPSHAcode 309 verification project (Hale et al., 2018). For a thorough comparison, we performed numerical computations using four dif- 310 ferent algorithms: Riemann sum, conventional Monte Carlo (MC), importance sampling (IS), and adaptive importance 311 sampling (AIS) using the VEGAS algorithm. These computations were conducted across various seismic source settings: 312 16 • www.bssaonline.org Volume 0 Number 0 Month Y ear Figure 2. Seimsic source geometry for the three comprehensive numerical examples to test our AIS PSHA framework. Figure 3. Benchmark PSHA curves for numerical examples 1 (areal source), 2 (linear fault source), and 3 (combined sources) 1)anarealsource(Area1),2)afaultsource(FaultA),and3)acombinationofanarealandtwofaultsources(Area1,FaultA, 313 FaultB). These examples’ geometry and seismic activity parameters are detailed in Figure 2 and Table 2. Figure 3 presents 314 thebenchmarkhazardcurvesfortheexamples.ThegroundmotionmodelfromSadighetal.(1997)wasemployed. 315 ForArealsource1,acircular-shapedarealsourcewitha100kmradiusisconsidered(Figure2).Theearthquakeactivity 316 rate of the source, 𝜈(𝑀>𝑚min), is 0.0395/year with 𝑚minand𝑚maxof 5.0 and 6.5, respectively. The earthquake occurrence 317 model is doubly-truncated exponential with Gutenberg-Richter 𝑏-value of 0.9. Seismogenic depth is 5 to 10 km from the 318 surface,andapointsourceisassumed. 319 Volume 0 Number 0 Month Y ear www.bssaonline.org •17 TABLE 2 Seismic source information for comprehensive numerical examples to test our AIS PSHA framework Area 1 Fault A Fault B Source Type Area Vertical Fault Vertical Fault Source Geometry Circle (R=100 km) Line (L = 50 km) Line (L = 85 km) Earthquake Occurrence Model Exponential Characteristic Characteristic b-value 0.9 0.9 0.9 Mmin 5.0 5.0 5.0 Mmax 6.5 6.75 7.0 M𝑐ℎ𝑎𝑟 - 6.5 6.75 𝜈(𝑀>𝑀min)(yr−1) 0.0395 - - slip rate (mm/yr) - 1 2 Seismogenic depth (km) 5-10 0-12 0-12 Rupture type Point Floating rectanglular rupture Floating rectanglular rupture FaultsAandBareverticalfaultsources(strike=90◦E,dip=90◦)withlengthsof50and85km,respectively.Theyextend 320 fromthesurfaceto12kmdepth.Theearthquakeactivitymodelsarechosentobecharacteristic(YoungsandCoppersmith, 321 1985)withthe 𝑏-valueof0.9,and 𝑚minissettobe5.0forbothfaultsources.ForFaultA,sliprate, 𝑀max,𝑀𝑐ℎ𝑎𝑟are1mm/year, 322 6.75,and6.5;andforfaultB,theyare2mm/year,7.0,and6.75.Infaultsources,rupturesareassumedtobefloatinginside 323 thefaultwiththerupturedimensionsfollowing 324 log10(𝐴)=𝑀−4 (35) log10(𝑊)=0.5𝑀−2.15 (36) log10(𝐿)=0.5𝑀−1.85 (37) ,where𝐴isruptureareainkm2,𝑊isrupturewidthinkm, 𝐿isrupturelengthinkm,and 𝑀isearthquakemagnitude.The 325 probabilityofexceedanceateighteengroundmotions(PGA),0.001,0.01,0.05,0.1,0.15,0.2,0.3,0.35,0.4,0.45,0.5,0.55,0.6, 326 0.7,0.8,0.9,and1.0g,areestimated. 327 For Riemann summation, we followed the calculation procedure specified in the PEER PSHA code verification project 328 (Hale et al., 2018). Thus, we used magnitude and source-to-site distance step sizes of 0.01 and 0.1 km, respectively. Note 329 that the number of grids of the source-to-site distance distribution ( ∼950) is significantly less than that of source location 330 probability, the uniformly distributed probability of event location across the seismic source region, with the spatial and 331 depthgridspacingof0.5kmand1km(𝜋×1002 0.5×0.5×6≈750,000)suggestedbyHaleetal.(2018).ThismakestheRiemannsum 332 computation ∼800timesmoreefficient.Forthegroundmotionrandomvariable( 𝜀),thegridstepsizewassetto0.01,with 333 minimumandmaximumvaluesof-6and6,respectively,beyondwhichthePSHAcurveshowsnegligiblechangeevenatan 334 annualexceedanceprobabilityof 10−8(BommerandAbrahamson,2006). 335 In the conventional MC approach, we generate a ground motion catalog based on the probability distributions of 𝑚,𝑟, 336 and𝜀.Fromthiscatalog,wecalculatetheannualexceedancefrequencybycountingtheinstanceswherethegroundmotion 337 18 • www.bssaonline.org Volume 0 Number 0 Month Y ear exceedsthespecifiedthreshold,asoutlinedinEq.(11).Thisfrequencyisthenconvertedtoanannualexceedanceprobability 338 usingEq.(6). 339 ForIS,weadoptuniformISdensitiesacrosstheentireintegrationrangeformagnitude,distance,and 𝜀followingJayaram 340 andBaker(2010).Thischoiceismadeduetotheabsenceofpriorinformationregardingwhichintegrationrangesignificantly 341 influencethehazardcalculation. 342 ForourAISPSHAwiththeVEGASalgorithm,wesettheinitialnumberofgridsfor 𝑚,𝑟,and𝜀to50(Eq.(28)).Thetotal 343 numberofsubgrids, 𝑁subgrid,ischosentobe10,000.Thelearningrate, 𝛼,isfixedat1.0.Algorithm1presentsthepseudocode. 344 Fromlines1to5,thegenerationindex 𝑡issettobezero,andthesamplingdensity 𝑞issettobeuniform.Themainalgorithm 345 loop is from lines 6 to 24. The loop is continued until there is no improvement in the coefficient of variation or there is no 346 previousgeneration(line6).Atline7, 𝑋𝑚,𝑋𝑟,𝑋𝜀aresampledfromthedistribution 𝑞(𝑡).Then,theprobability ̂𝜆isestimated 347 in lines 8-10, and its coefficient of variation is also calculated in line 11. From lines 12 to 23 the sampling function 𝑞(𝑡)is 348 updated.Thestepsarerepeatedover 𝑚,𝑟,and𝜀(lines12-15).Fromlines16to20,thecontributionofeachgridiscalculated, 349 atline21,itissmoothedanddampened,and 𝑞(𝑡)isupdatedbysubdivisiondependingon 𝑑𝑖andrestorethenumberofgridsto 350 theoriginalnumber, 𝑁,inline22.Updated 𝑞(𝑡+1)isobtainedbymultiplying 𝑞(𝑡+1) 𝑚,𝑞(𝑡+1) 𝑟,and𝑞(𝑡+1) 𝜀inline23,andthewhile 351 loopisendedbyincreasethegenerationindex 𝑡.Whenmainalgorithmloopisterminated,itreturnsthehazardestimate ̂𝜆 352 andproposedoptimaldensity 𝑞inline25. 353 We assessed the accuracy of conventional MC, IS, and AIS probability estimates through the standard deviation of the 354 relativeerrorwithrespecttothebenchmarkcurve,calculatedusingthefollowingformula: 355 𝜎=√ √√√√√1 𝑁𝑁∑ 𝑛=1𝑒2 𝑛−⎛ ⎜ ⎝1 𝑁𝑁∑ 𝑛=1𝑒𝑛⎞ ⎟ ⎠2 , where𝑒𝑛is the relative error of the 𝑛th MC estimate, and 𝑁is the total number of MC exceedance probability estimates. 356 Since MC, IS MC, and VEGAS AIS estimates are unbiased (Eq. (12), Eq. (23)), the sum of 𝑒𝑛when𝑁→∞is theoretically 357 zero.Therelativeerror, 𝑒𝑛,isdefinedas: 358 𝑒𝑛=𝑝𝑛−𝑝𝑟𝑒𝑓 𝑝𝑟𝑒𝑓 ,where𝑝𝑛istheexceedanceprobabilityestimatedat 𝑛thMCestimate,and 𝑝𝑟𝑒𝑓isthebenchmarkprobability. 359 TheanalyseswereconductedinaPython3.11onanIntelCorei7-137002100MHzprocessorwith64GBRAM. 360 Example 1 : Areal source 361 Weconsiderthecircularareasourcewiththesiteatthecircle’scenterinFigure2.Forthisexample,weshowdifferentcompu- 362 tationalperformancestoconductPSHAinFigure4.Ournumericalexperimentsshowthatcomputationaltimesandstandard 363 Volume 0 Number 0 Month Y ear www.bssaonline.org •19 Algorithm 1 VEGASadaptiveimportancesamplingPSHApseudocode [Parameters and Functions] 𝑎 Groundmotionintensityofinterest(e.g., 𝑎=0.1g) 𝑁 Numberofgrids(e.g., 𝑁=50) 𝑁subgridNumberofsub-grids(e.g., 𝑁subgrid=10,000) 𝑁𝑠NumberofMCsamples(e.g., 𝑁𝑠=1,000) ⃗𝑋𝑚Magnitudesamplevector( ⃗𝑋𝑚,𝑖=𝑖thelementof ⃗𝑋𝑚) ⃗𝑋𝑟Distancesamplevector ⃗𝑋𝜀Groundmotionrandomvariablesamplevector ⃗𝑋 [⃗𝑋𝑚;⃗𝑋𝑟;⃗𝑋𝜀] 𝐺(⋅)GroundMotionModel 𝐼(⋅)Indicatorfunction 𝑓𝑋(⋅)Originalsamplingdistribution 𝑞(𝑡)(⋅)Proposedsamplingdistributionat 𝑡thiterationstep ̂𝜆(𝑡)Hazardestimateat 𝑡thiterationstep COV(𝑡)COVofthehazardestimateat 𝑡thiterationstep 𝜖 pre-definediterationstoppingcriteria(e.g.COVof0.02) 𝑖 Indexforpartitionedgrids( 𝑖=1,2,⋯,𝑁) 𝑗 Indexforsamples( 𝑗=1,2,⋯,𝑁𝑠) [Algorithm] 1:𝑡=0 2:𝑞(𝑡) 𝑚=U([𝑚min,𝑚max]) 3:𝑞(𝑡) 𝑟=U([𝑟min,𝑟max]) 4:𝑞(𝑡) 𝜀=U([𝜀min,𝜀max]) 5:𝑞(𝑡)=𝑞(𝑡) 𝑚𝑞(𝑡) 𝑟𝑞(𝑡) 𝜀 6: whileCOV(𝑡−1)<𝜖ort<2: 7: ⃗𝑋={𝑋𝑚,𝑋𝑟,𝑋𝜀}𝑁𝑠,(𝑋𝑚∼𝑞(𝑡) 𝑚, 𝑋𝑟∼𝑞𝑟, 𝑋𝜀∼𝑞(𝑡) 𝜀) 8: ⃗𝐻(𝑡)={𝐼(𝐺(𝑥)>𝑎|𝑥)𝑓𝑋(𝑥)|𝑥∈⃗𝑋} 9: ⃗𝜆(𝑡)=⃗𝐻(𝑡)∕𝑞(𝑡)(⃗𝑋) 10: ̂𝜆(𝑡)=∑𝑁𝑠 1𝜆(𝑡) 𝑖∕𝑁𝑠 11: COV(𝑡)←Eq.(25) 12: for𝑢in[𝑚,𝑟,𝜀]: 13: if𝑢=𝑚:(𝑣,𝑤)←(𝑟,𝜀) 14: if𝑢=𝑟:(𝑣,𝑤)←(𝜀,𝑚) 15: if𝑢=𝜀:(𝑣,𝑤)←(𝑚,𝑟) 16: for𝑖in{1,2,...,𝑁}: 17: for𝑗in{1,2,...,𝑁𝑠}: 18: if𝑢𝑖−1<𝑋𝑢,𝑗<𝑢𝑖:𝐻𝑖←𝐻𝑖+(𝐻2 𝑗∕(𝑞𝑣(𝑋𝑣,𝑗)⋅𝑞𝑤(𝑋𝑤,𝑗)) 19: 𝐻𝑖←√ 𝐻𝑖 20: 𝑑𝑖←(𝐻𝑖∆𝑥𝑖∑ 𝑖𝐻𝑖∆𝑥𝑖) 21: 𝑑𝑖←smoothed,dampened 𝑑𝑖(Eq.(33),(34)) 22: 𝑞(𝑡+1) 𝑢 ←subdivisionandrestoration(Eq.(29)-(31)) 23: 𝑞(𝑡+1)←𝑞(𝑡+1) 𝑚𝑞(𝑡+1) 𝑟𝑞(𝑡+1) 𝜀 24: 𝑡←𝑡+1 25: return̂𝜆(𝑡−1),𝑞(𝑡−1) 20 • www.bssaonline.org Volume 0 Number 0 Month Y ear deviations(“accuracy”)arelinearlycorrelatedinlogarithmicscale,inagreementwiththetheory,becausetherequirednum- 364 berofMCsamples(linearlyproportionaltothecomputationaltime)isinverselyproportionaltothesquareofthestandard 365 deviation(Eq.(16)). 366 While the computation time and standard deviation vary depending on the target ground motion of interest, our results 367 show our AIS PSHA generally outperforms the other numerical techniques. AIS PSHA becomes increasingly efficient for 368 larger ground motions. At a low target ground motion of 0.05g, the computation time to achieve a 2 % standard deviation 369 is 0.01 seconds for AIS, while it takes 0.05 and 0.09 seconds for conventional MC and IS estimates, respectively, indicating 370 AIS is 4.8 and 8.5 times faster. However, computational efficiency becomes extreme in high-ground motions. At 1.0 g, the 371 computationtimetoachieve2%standarddeviationisestimatedtobe0.02secondsforAIS,whileittakes127and1.4seconds 372 forconventionalMCandISestimates.Inthiscase,AISis7,800and70timesfasterthanconventionalMCandIS,respectively. 373 Also,notethatfor2%standarddeviationcase,AISis >105fasterthanRiemannsum.WealsonotethatAISoutperformsIS 374 inallgroundmotionrangesbyafactorof8to70toachievea2%standarddeviation.However,thisisnotalwaysguaranteed 375 because AIS takes 𝑡times more computation time than IS with the same 𝑁due to the 𝑡iterations to find the optimal IS 376 density.Thus,ourfindingsimplyAISPSHAwiththeVEGASalgorithmcanfind(closeto)optimalISdensityquickly. 377 Itisalsonoteworthythatwithfixednumberofsamples(i.e.,constantcomputationalburden),AIS’sestimationaccuracyis 378 quitesimilaracrossdifferenttargetgroundmotions,whileconventionalMC’saccuracydecreasessharplyforhigherground 379 motions(Figure5).Infact,ifthehazardcurveisexponential(MarzocchiandJordan,2017),wecanshowconventionalMC 380 decreasesitsefficiencyalsoexponentially(Eq.(15)).Forexample,theerrorinthehazardcurvefor1.0gwillincreaseto100% 381 ifitinitiallywas1%for0.001gwithexceedancefrequencyof0.9/year(when 𝜈=1/year).Fromournumericalexperiments, 382 we found errors grow from 0 to 5.60 % with 𝑁convMC=107MC samples for these ground motions. In contrast, for a similar 383 computationtime,wefoundAIShaderrorsrangingfrom0.13to0.64%,forthisquitedifferentgroundmotionlevels(Figure 384 5). 385 ThisfindingiskeyforPSHAasthecomputationalbottleneckisatthehighestgroundmotionintensity.Usingconventional 386 MC PSHA, the hazard analyst has no choice but to largely increase the number of samples to estimate hazard accurately 387 at high ground motions even though such a large number would not be necessary for low ground motions. This makes 388 theconventionalMChighlyinefficient,andconsequently,theefficiencyoftheconventionalMCforlowergroundmotions 389 cannotbeconsideredarealadvantageforPSHA.AISPSHAovercomesthisproblembyadoptingdifferentoptimaldensities 390 atdifferentgroundmotionintensities,makingthecomputationalburdentoachievethesimilarestimationerroralmostflat 391 foranygroundmotionintensity,asshowninFigureS1.IntermsoftheaccuracyofAISPSHA,thePEERPSHAverification 392 projectsuggestsastrictacceptableerrorrangeof5%forreliablePSHAcomputationcodes.AISPSHAachievesitwithonly 393 N∼10,000pergroundmotion(FigureS1). 394 Volume 0 Number 0 Month Y ear www.bssaonline.org •21 AnotherkeyadvantageofourproposedAISPSHAistheco-productionofdeaggregationcurvesatnoextracomputational 395 cost.FromEq.(26),optimalISdensitiesaretheoreticallyequivalenttohazarddeaggregationdistributions.WeshowAIScan 396 find close-to-optimal IS densities and thus closely match hazard deaggregation curves. We found the benchmark marginal 397 distributions of hazard deaggregation obtained from Riemann sum closely match the iterated IS density from AIS PSHA 398 (Figure6).WeusedtheKolmogorov-Smirnov(K-S) 𝐷statistic(Kolmogorov,1933)toquantifytheirsimilarities.K-S 𝐷statistic 399 measuresthemaximumdifferencebetweentwocumulativedistributionfunctions(CDF).IftwoCDFareidentical, 𝐷iszero, 400 anditsmaximumpossiblevalueisone. 𝐷closetozeroindicatesthattwoprobabilitydistributionsaresimilar.Wecalculated 401 𝐷for𝑚,𝑟, and𝜀at all the ground motion intensities (Figure S2). We found the maximum 𝐷values for𝑚,𝑟,𝜀were 0.032, 402 0.113,and0.092,respectively,andtheminimumvalueswere0.019,0.026,and0.017,indicatingastrongresemblancebetween 403 thetwodistributions. 404 Wealsoestimatedthedifferencesofmeanvaluesfromproposed 𝑞∗andRiemannsumhazarddeaggregation(FigureS2). 405 NotethatthedeaggregationdistributionsandoptimalISdensitiesvaryfordifferentgroundmotionlevels.Thus,theirmean 406 values also vary. We found the maximum relative differences were 2.5 %, 22.6 %, and 18.5 % for 𝑚,𝑟,𝜀, and the minimum 407 differences were 1.7 %, 1.0 %, and 4.3 %, respectively. Also, the maximum absolute differences were 0.14, 3.9 km, and 0.17 408 for𝑚,𝑟,𝜀, and the minimum differences were 0.09, 0.7 km, and 0.0009. The maximum relative difference in distance ( 𝑟) 409 appears for deaggregation distributions at 0.25 g, where the mean distance obtained from hazard deaggregation is 14.3 km 410 andthatfrom 𝑞∗is17.5km.Giventhatwearetypicallyinterestedindistancerangesontheorderoftensofkilometers(e.g., 411 0-15 km, 15-25 km, 25-50 km, etc.) rather than a single value (U. S. Nuclear Regulatory Commission, 2007), this difference 412 isnotcrucialindeterminingthecontrollingearthquakeforcriticalinfrastructures.Inaddition,ISdensitiesstillshowsmall 413 K-S𝐷statisticsandtheirmodealmostmatcheseachother,asshowninFigure6. 414 Example 2 : Fault source 415 Weconsidera50km-lengthverticalfault25kmawayfromthesite(Figure2).Weadoptedafinite-dimensionrupturemodel, 416 whichresultsinadistancedistributiondependentonmagnitude.Wefoundthisdependencycandiminishtheperformance 417 of VEGAS AIS because the VEGAS algorithm assumes the independently distributed random variables. For example, for a 418 ground motion intensity of 0.05 g, AIS is slower than conventional MC (Figure 7), while AIS outperformed conventional 419 MC at the same ground motion intensity when point source assumption was made, as seen earlier (Figure 4). However, as 420 thegroundmotionintensityincreases,thecomputationalgapbetweenthetwomethodsbecomessmallerrapidlyandcloses 421 at 0.2g. For higher ground motions, AIS outperforms conventional MC. At 1.0 g, for a 5% standard deviation, AIS, IS, and 422 conventionalMCtake0.08,3.7,and166seconds,respectively,i.e.,AISMCis48and2,162timesfaster. 423 We also present the accuracy of AIS PSHA in comparison to the true hazard calculated by Riemann summation (Figure 424 S3). We observed AIS PSHA estimates approximate the true hazard curve within an acceptable range (5 % error) when 𝑁 425 22 • www.bssaonline.org Volume 0 Number 0 Month Y ear Figure 4. Areal Source Example: Standard deviation of the conventional MC (red), IS MC (blue), AIS MC (black) estimates as functions of computation time at ground motion of 0.05 g, 0.2 g, 0.5 g, and 1.0 g. The computation time required for Riemann summation is presented as green vertical dashed line. At low ground motion, conventional MC outperforms IS and AIS MC, however, as the target ground motion increases, the performance of AIS make a dramatic improvement in terms of both accuracy and computational cost. is greater than ∼50,000. Also, note that though we assumed the independently distributed optimal density, AIS PSHA still 426 gives an unbiased hazard estimate that will converge to the true hazard with a sufficiently large number of samples due to 427 thenatureoftheimportancesampling(FigureS3(d)). 428 WealsocomparedhazarddeaggregationandtheiteratedISdensity(Figure8)andshowedtheycloselymatcheachother 429 eventhoughthedistancedistributiondependsonthemagnitudeinthiscase.ItisalsonoteworthytheiteratedISdensitycan 430 evenreproducecomplexdensitieswithdiscontinuitieslikethelargejumpwithinmagnitudedistribution(for 𝑚=6.25)dueto 431 theuseofacharacteristicearthquakeoccurrencemodel.TheK-SDstatisticandmeandifferenceofthetwodistributionsare 432 alsopresentedinFigureS4.ThemaximumvaluesofDin 𝑚,𝑟,and𝜀are0.30,0.68,and0.13,respectively,andtheminimum 433 valuesare0.02,0.31,and0.04.Wefoundthatthelargestdiscrepanciesoccurinthemagnitudedistribution,buterrorscanbe 434 considered negligible as the mean magnitude difference is within 6 % error. We also note considerable discrepancies in the 435 distancedistributionshape(seeK-SDstatisticsofFigureS4)ascurveswithconcentratedprobabilitiesinnarrowrangesare 436 Volume 0 Number 0 Month Y ear www.bssaonline.org •23 Figure 5. Standard deviation of MC estimates with the ground motion for conventional (red) MC ( 𝑁convMC=107) and AIS (black) with the similar computation time. The standard deviation of error exponentially increase with the ground motion in conventional MC, while that of AIS remain constant. hardertoestimateforAIS(Figure8).However,themeandistanceisstillwithina5%errorrange.The 𝜀isgenerallyingood 437 agreementacrossallgroundmotionintensities.Weobservehigherdifferencesatlowergroundmotionbecausethemean 𝜀 438 is close to zero. For example, the calculated mean 𝜀at ground motion intensity of 0.01 g is 0.022 and 0.007, which is not a 439 largedifferenceinpractice. 440 Themaximumrelativedifferencesinmeanvalueswerefoundtobe6.2%,5.3%,and129%for 𝑚,𝑟,𝜀,andtheminimum 441 differences were 0.006 %, 3.9 %, and 2.2 %, respectively. Note that the 129 % of 𝜀case is corresponding the case where the 442 mean𝜀isclosetozero.Therelativedifferenceappearstobeslightlyhigherthanpreviousarealsourceexample,however,the 443 absolutedifferenceisstillremaintobesignificantlysmall,themaximumabsolutedifferenceswere0.36,1.3km,and0.20for 444 𝑚,𝑟,𝜀,andtheminimumdifferenceswere0.0004,1.0km,and0.07. 445 Example 3 : Combined sources 446 InPSHA,weoftenhavemultipleseismicsources.Weconsideronearealandtwofaultsourcesaroundthesitetorepresentthis 447 case.Thisapplicationpositsadifferentmathematicalproblemthantheprevioustwoexamplesbecausewemustintroduce 448 anadditionalvariabletoformulateAIS. 449 First,theprobabilityofearthquakeoccurrenceat 𝑖thseismicsourcecanbedefinedas: 450 𝑃(𝑆=𝑖)=𝜈𝑖 ∑𝑁𝑠 𝑖=𝑗𝜈𝑗 ,where𝑁𝑠isthenumberofseismicsources( 𝑁𝑠=3inthisexample),and 𝜈𝑖istheannualearthquakeoccurrencerateof 𝑖th 451 seismicsource.BecausethediscreterandomvariablescannotbeusedinAIS,wedefineacontinuousrandomvariableand 452 itscorrespondingprobabilitydensityfunctionas 453 24 • www.bssaonline.org Volume 0 Number 0 Month Y ear Figure 6. Areal Source Example: The convergence of 𝑚,𝑟, and𝜀iterated IS densities density (red) derived in AIS algorithm (N = 100,000) to marginal distributions of hazard deaggregation (black) at ground motion intensities of 0.001 g, 0.01 g, and 0.5 g. 𝑓𝑆(𝑠)=⎧ ⎪⎪ ⎨⎪⎪ ⎩0𝑠<0 𝑃(𝑆=𝑖)if𝑖−1≤𝑠 < 𝑖(𝑖=1,2,...,𝑁𝑠), 0𝑁𝑠≤𝑠(38) 𝑓𝑆(𝑠)is a piece-wise constant function where the heights are proportional to the corresponding sources’ earthquake 454 occurrencerates.Weintroduce 𝑓𝑆(𝑠)intothePSHAintegrationandobtain 455 Volume 0 Number 0 Month Y ear www.bssaonline.org •25 Figure 7. Fault Source Example: Standard deviation of the conventional MC (red), IS MC (blue), AIS MC (black) estimates as functions of computation time at ground motion of 0.05 g, 0.2 g, 0.5 g, and 1.0 g. Λ(𝑋>𝑎)=⎛ ⎜ ⎝𝑁𝑠𝑟𝑐∑ 𝑖𝜈𝑖⎞ ⎟ ⎠∫ 𝑠∫ 𝜀∫ 𝑅∫ 𝑀𝐼(𝑥>𝑎|𝑠,𝑚,𝑟,𝜀)𝑓𝑀,𝑅,ℰ(𝑚,𝑟,𝜀|𝑠)𝑓𝑆(𝑠)𝑑𝑚𝑑𝑟𝑑𝜀𝑑𝑠 (39) ,whichistheintegralversionof Eq.(1)andEq.(5).SolvingEq.(39)positscomputationalchallengesthanthesinglesource 456 problemsbecause 𝑓𝑆(𝑠)havethelargejumpsat 𝑖=1,2,…,𝑁𝑠−1.Inaddition,weintroduceadditionaldependenciesinthe 457 multi-source case because 𝑚and𝑟depend on 𝑠, further diminishing the VEGAS algorithm’s effectiveness. Thus, we test 458 twoAISPSHAapproachesforthiscase:1)fullAISapproachutilizingtheEq.(39)and2)partialAISapproach,whichisthe 459 simplesummationofsingle-sourceAISPSHAcurves. 460 WecomparedthecomputationalperformanceofMC,IS,andfullandpartialAIS(Figure9).Likethepreviousexamples, 461 conventional MC is faster for low ground motions, e.g., 0.05 g, but not high ones. For example, to achieve the 𝜎of 2 % at a 462 ground motion of 0.8 g, conventional MC, IS, and full and partial AIS take 73.8, 10.0, 0.77, and 0.13 seconds, respectively, 463 i.e.,partialAISisthemostefficientalgorithm, ∼583,79,and6timesfasterthanconventionalMC,uniformIS,andfullAIS, 464 respectively. 465 26 • www.bssaonline.org Volume 0 Number 0 Month Y ear Figure 8. Fault Source Example: The convergence of 𝑚,𝑟, and𝜀iterated IS densities density (red) derived in AIS algorithm (N = 1,000,000) to marginal distributions of hazard deaggregation (black) at ground motion intensities of 0.1 g, 0.5 g, and 0.8 g. As stated earlier, full AIS is less efficient than partial AIS because of the additional dependencies and jumps introduced 466 by𝑓𝑆(𝑠). Figures S5 and S6 present the accuracy of both AIS approaches with different sample sizes. The contrast between 467 Figures S5 (c) S6 (c) shows that partial AIS has smaller error than full AIS even with fewer samples (150,000 < 500,000). 468 Wedon’tshowthecomparisonforthepartialAISongroundmotionsgreaterthan0.8gbecausefaultBreachesanumerical 469 instabilityduetoitslowexceedanceprobability(< 10−12/yr)(FigureS5). 470 Volume 0 Number 0 Month Y ear www.bssaonline.org •27 Figure 9. Combined Sources Example: Standard deviation of the conventional MC (red), IS MC (blue), full AIS MC (black), and partial AIS MC (gray) estimates as functions of computation time at ground motion of 0.05 g, 0.25 g, 0.4 g, and 0.8 g CONCLUSION 471 We proposed a novel computational framework for PSHA based on an implementation of the VEGAS algorithm. Through 472 comprehensive testbeds, we investigated the computational performance of this new approach, covering widely adopted 473 source types in PSHA practice: 1) areal source, 2) fault source, and 3) combined sources, including both areal and fault 474 sources. 475 Wecomparedtheproposedmethodtothreeexistingcomputationalframeworks:a)Riemannsum,whichhasexponentially 476 increasing computational costs for finer grid sizes; b) conventional MC, which requires a substantially long catalog with 477 synthetic earthquakes, particularly for large seismic intensities; and c) importance sampling (IS) with simple (uniform) IS 478 distributions. 479 OurfindingsindicatethatAISPSHAoutperformsallothercomputationalframeworks.AISPSHAcandramaticallyreduce 480 computationaltimesbyfactorsupto >105comparedtotraditionalRiemannsummation.AISPSHAwasalso 103fasterthan 481 28 • www.bssaonline.org Volume 0 Number 0 Month Y ear theconventionalMCwhilemaintaininga2%standarddeviationoferror.Additionally,AISPSHAwasupto70timesfaster 482 thanISPSHA,demonstratingthattheVEGASalgorithmcanapproximateoptimalISdistributionsquicklyandwell. 483 We showed that AIS PSHA requires a similar computation time for any ground motion, making its application to larger 484 ground motions with low probability substantially more efficient than conventional MC. In PSHA practice, the computa- 485 tional demands are dominated by large ground motions that need substantially more samples than lower ground motions. 486 Incontrast,AISPSHAonlyrequiresasimilarlylownumberofsamplesforallgroundmotionlevelsbecauseitfindsoptimized 487 ISdistributionsforeach. 488 WealsoshowedthatAISPSHAfindsapproximateddeaggregationcurvesatnoextracomputationalcostbasedontheoret- 489 ical insights showing that optimal IS densities are equivalent to deaggregation distributions. We show empirically that the 490 hazard deaggregation and iterated IS densities from AIS PSHA are fairly similar by comparing the statistical properties of 491 thetwodistributions,e.g.K-SDstatistics <0.113andmeanvaluesdifferencesof <4.3%. 492 Our study indicated that our AIS PSHA implementation works extremely well for point sources, where the magnitude 493 anddistanceareindependentrandomvariables.Wealsoappliedthealgorithmtoproblememployingfinite-rupturemodel, 494 and it still outperformed the pre-existing algorithms up to by factors of >2,000. It appears to be less effective than point 495 sourceexamplebecauseitintroducesadependencybetweenthemagnitudeanddistance.OurimplementationoftheVEGAS 496 algorithmusesindependentrandomvariablesfortheISdistribution,thus,makingitlesseffectiveforthefinite-rupturecase. 497 However, future implementations can also consider other versions of AIS algorithms where the variables are correlated; 498 thoughcomputationaldemandsfortheAISiterationwouldincrease. 499 Inthecaseofcombinedseismicsources,weproposedtwostrategies:1)incorporatingthesourcerandomvariableintothe 500 AIS PSHA framework and 2) simple summation of AIS PSHA curves for individual sources. Both strategies outperformed 501 the traditional methods up to by a factor of ∼580. However, the second strategy was more efficient than the first one by a 502 factorof6.Thefirststrategyaddeddependenciestotheseismichazarddistribution,makingtheVEGASlesseffective. 503 In sum, AIS can be applied to any PSHA computation, leading us to expect widespread application of the method. 504 Specifically,weconsiderourproposedAISPSHAtobesignificantlybeneficialforlarge-scaleprojectsthatinvolvenumerous 505 logictreebranchesandhaveextremecomputationaldemands. 506 DATA AND RESOURCES 507 The source code for computing PSHA using the framework explained in this paper is available at 508 https://github.com/sehoung/ais_psha. Figures were created using Matplotlib (Hunter, 2007) and Microsoft PowerPoint 509 (http://office.microsoft.com;lastaccessedJune2024). 510 Volume 0 Number 0 Month Y ear www.bssaonline.org •29 DECLARATION OF COMPETING INTERESTS 511 Theauthorshavefiledapatentapplicationrelatedtothemethoddiscussedinthispaper.Thisapplicationisrelevanttothe 512 researchpresentedandcouldrepresentapotentialfinancialinterest. 513 ACKNOWLEDGMENTS 514 WeacknowledgethefinancialsupportbytheTandonSchoolofEngineeringattheNewYorkUniversityandthedepartmentofCiviland 515 EnvironmentalEngineeringattheUniversityofCalifornia,Berkeley. 516