Losch, Martin

Schröter, Jens

Sloyan, B.

Yaremchuk, A.

Yaremchuk, M.

English

Electronic Publication Information Center

On the Estimation of Posterior Error BarsMartin Losch1, Alexey Yaremchuk1,2, Maxim Yaremchuk3, Jens Schröter1, Bernadette Sloyan11Alfred-Wegener-Institute for Polar- and Marine Resarch, Bremerhaven, Germany2Andreyev Acoustics Institute, Moscow, Russia3International Pacific Research Center (IPRC), University of Hawaii, USA(email: mlosch@awi-bremerhaven.de)IntroductionPosterior error estimates are important parts of an inverse method solution. Unfor-tunately, their evaluation can become computationally quite demanding since thistask always involves some sort of inverse of the Hessian of the problem. However,some observables for which the error estimates are required can be represented bya relatively small part of the Hessian spectrum and thus do not need its full inversefor estimating their errors.We propose a polynomial approximation of the Hessian spectrum which enablesus to obtain reliable error estimates for some observables. We compare the resultsto “standard” methods and show examples of applications for large problems.Methods1. Full inversionThe inverse of the Hessian of the cost function is the covariance of the deviation ofthe control variables ϕ from the optimal state ϕ∗. An observable Φα, α=1,2,... canbe expanded into powers of (ϕ−ϕ∗):One finds that near the optimum state the covariance of the observables can beexpressed as φαTH-1φβ .2. Conjugate gradientsWhen the dimension D of the problem is very large, it is still possible to calculatethe vector Hψ for any ψ. In order to find φαTH-1φβ, the equation HXβ = φβ issolved for Xβ with a conjugate gradients method to obtain the covariance φαTXβ.3. Polynomial ApproximationThe covariance φαTH-1φβ is written in a more symmetric form ψαTψβ withψ = H-1/2φ. The latter can be expressed as a series of Chebyshev's PolynomialsUn(x):with un = Un (1 - 2H) φ. Truncuation of the series corresponds to smoothing thespectrum of the Hessian H. Let εk = sin2(θk/2) be the eigenvalues of H and ek theeigenvectors, then the spectral functionis the spectral function of the Hessian δ(θ−θk) weighted by (φαTek)(φβTek). Truncu-ation results again in smoothing the spectrum. From this spectrum we can decidewhat part of the spectrum of H is needed for the correct evaluation of φαTH-1φβand choose the inverting window accordingly.4. ComparisionFor the Southern Ocean WOCE section SR3 a small geostrophic toy model with136 control variables (tracers and bottom velocities) illustrates the three methodsdescribed. Figures 1 and 2 show the temperature and velocity field after inversion.the posterior errors and their correlation of the total mass and heat transportsthrough the section obtained by the three different methods show a reasonableagreement (Table 1)Table 1:Comparision of errors obtained by the three different methodsMethod # of iterations Mass TransportError [Sv]Heat TransportError [PW]Correlation ofErrorsFull inversion           32          1.6          0.27conjugategradients       520           32          1.7          0.17polynomialapproximation       500           32          1.4          0.28Φα Φα ϕ∗( ) φαT ϕ ϕ∗–( ) …+ +=H12--– φ 16pi-----n 1+4 n 1+( )2 1–------------------------------- un   n 0=∞∑=ραβ δ θ θk–( ) φαTek( ) φβTek( )k 1=D∑=θ2pi-- nθ( ) φαTun 1–β( )( )sink 1=∞∑sin=10−7 10−6 10−5 10−4 10−3 10−2 10−1 1000102030405060708090100110EigenValues (Normalized)Number of EigenValues (in %)0 20 40 60 80 100 12010−710−610−510−410−310−210−1100EigenValues (Normalized)Number of EigenValuesFigure 1: Temperature of toy model for SR3. Note the coarse grid.Figure 4: Eigenvaluespectrum of the Hessian matrix of the toymodel. The condition number is 3x10-6Figure 5: Same eigenvaluespectrum of the Hessian matrix of thetoy model for SR3 as in Figure 4. The representation has beenchanged to match that of figures 6-8 (lower panels). On theabscissa, there are the normalized magnitude of the eigenvalues, onthe ordinate the index number of the eigenvalues relative to thedimension of the problem (D = 136)Figure 2: Velocity normal to section (cm/s) of the toy model of SR3.Figure 3:AWI10−5 10−4 10−3 10−2 10−1 100−50−40−30−20−10010EigenValues (Normalized)Magnitude (dB) Spectral Density Relative to Uniform Fourier DistributionEval_funct_U: 05/26/99 08:47:38, /tphs1/user5/mlosch/conferences/Quebec/case2/isqrSpd500.datChebyshev (50 dB) smoothing, Nspd=500, Nvac=2, Mplt=350makeISQ −> Emin:  1.903e−05, Nfilter: 500, acc: 0.0110−5 10−4 10−3 10−2 10−1 1000102030405060708090100110Cumulative Spectral Distribution Versus EigenvaluesDistribution function (percent)EigenValues (Normalized)−15 −10 −5 0 5 10 15 2012345Water MassMass Flux (× 106 m3 s−1)Indian 32° SSurface/Thermocline WaterIntermediate WaterUpper Deep WaterLower Deep WaterBottom Water10−7 10−6 10−5 10−4 10−3 10−2 10−1 100−50−40−30−20−10010EigenValues (Normalized)Magnitude (dB) Density of Eigenvalues Relative to Uniform Fourier DistributionHessian_spd_U: 05/08/99 00:10:32/home/awi1/mlosch/FLAME/24N/24N_03−Jan−1911/GRACE10/HeSp4096.datChebyshev (40 dB) smoothing, Nspd=4096, Nvac=2, Mplt=35010−7 10−6 10−5 10−4 10−3 10−2 10−1 1000102030405060708090100110Cumulative Spectral Distribution Versus EigenvaluesDistribution function (percent)EigenValues (Normalized)10−7 10−6 10−5 10−4 10−3 10−2 10−1 100−50−40−30−20−10010EigenValues (Normalized)Magnitude (dB) Spectral Density Relative to Uniform Fourier DistributionVac_spd_U: 05/26/99 08:33:38, /tphs1/user5/mlosch/conferences/Quebec/case2/VacSpd4096.datChebyshev (50 dB) smoothing, Nspd=4096, Nvac=2, Mplt=350makeISQ −> Emin:  1.903e−05, Nfilter: 500, acc: 0.0110−7 10−6 10−5 10−4 10−3 10−2 10−1 1000102030405060708090100110Cumulative Spectral Distribution Versus EigenvaluesDistribution function (percent)10−7 10−6 10−5 10−4 10−3 10−2 10−1 100−50−40−30−20−10010EigenValues (Normalized)Magnitude (dB) Density of Eigenvalues Relative to Uniform Fourier DistributionHessian_spd_U: 05/26/99 08:31:51/tphs1/user5/mlosch/conferences/Quebec/case2/HeSp4096.datChebyshev (40 dB) smoothing, Nspd=4096, Nvac=3, Mplt=35010−7 10−6 10−5 10−4 10−3 10−2 10−1 1000102030405060708090100110Cumulative Spectral Distribution Versus EigenvaluesDistribution function (percent)EigenValues (Normalized)Figure 6:• Upper panel: smoothed spectral density of Hessian of toy model for SR3.Number of coefficients un(x) = 4096.• Lower panel: cumulative distribution of the same Hessian spectrum.Going from right to left you can see how many eigenvalues (in %) areresolved. Some of the peaks in the upper part of the spectrum are so sharpthat they cannot be plotted, but their contribution shows in the cumlativedistribution.Figure 7:• Upper panel: smoothed weighted spectral density of the toy model formass (cyan line) and heat transport (purple filled contour). Number ofcoefficients un(x) = 4096.• Lower panel: cumulative distribution of the weighted spectrum. Since thedistribution function drops to zero at about 2x10-5 of the maximum eigen-value for both properties (cyan line: volume transport, purple filled con-tour: heat transport), the inverting window (red line) can be chosen suchthat only 500 coefficients un(x) are needed for the inversion, thus reducingthe computational cost considerably. As in figure 6 some eigenvalues inthe upper part of the spectrum cannot be seen in the upper panel but theyshow in the cumulative distribution.Figure 8Same weighted spectral density (upper panel) and cumulative spectral distri-bution (lower panel) for mass (cyan line) and heat transport (purple filledcontour) as in figure 7, but with a smaller number of coefficients un(x) =500. The smaller number of coefficients leads to a higher degree of smooth-ing of the spectral density.1. Box model of the Southern OceanAn inverse box model of the Southern Ocean estimates heat and volume fluxes perlayer through hydrographic sections. The control space has 1600 dimensions. Thismakes the full inversion of the Hessian possible. As an example of the results, themass flux and its errors for the zonal Indian Ocean section along 32oS (Ind32) areshown in figure 9. Table 2 lists the correlation of the flux errors of the differentwater massesTable 2: Correlation of mass flux errors for Indian 32oS (SW = Surface/Ther-mocline Water, IW = Intermediate Water, U/LDW = Upper/Lower Deep Water,BW = Bottom Water)SW IW UDW LDW BWSW     1     0.07     0.06   -0.05   -0.21IW     0.07     1   -0.51   -0.09     0.05UDW     0.06   -0.51     1   -0.47   -0.07LDW   -0.05   -0.09   -0.47     1   -0.59BW   -0.21     0.05   -0.07   -0.59     1Figure 9:Mass Transports (solid) and posterior errors (dash-dotted), divided intowater masses for a zonal section in the Indian Ocean at 32oS.Examples of Applications2. Nonlinear section inverse of SR3A nonlinear geostrophic inverse model similar to the one used as the toy model(see above) is applied to the meridional Southern Ocean WOCE section SR3between Tasmania and the Antartic continent, to estimate mass and propertytransports. As the number of control variables is quite large (~17000), the Hessiancannot be calculated nor inverted explicitly. The errors  in Table 3 are estimatedusing the conjugate gradients approach, where for each variance estimate severalthousand iterations were necessary for the algorithm to converge.Table 3: Transport estimates through SR3Mass [106 m3/s]Heat [PW]Salt[108 m3/s psu]Oxygen[Mmol/s]Silicate[102 kmol/s]   146 +/- 29  1.72 +/- 0.33      50 +/- 10      34 +/- 7      82+/-253. Nonlinear section inverse -- a twin experimentThe same model as in the previous section but augmented with satellite altime-try data is used to invert an artifical hydrographic section data set produced by ageneral circulation model (1/3o North Atlantic Model of the FLAME group).The zonal section is located in the North Atlantic along 24.5oN. Again the num-ber of control variables is fairly large (~12000) and the Hessian cannot be com-puted completely. Here for the variance estimation of the mass and heattransports, the polynomial approximation approach is employed. Figure 10 and11 show the spectral density of the Hessian of the model and the weighted spec-tral density for mass and heat transports through one small part of the section(Florida Strait). Note that the number of coefficients used to represent the Hes-sian spectrum is not bigger than for the toy model, although the problem ismuch larger, neither does the number of coefficients needed for the inversiondepend directly on the dimension of the problem. The transports and errorsobtained by this method are for mass 28.2 +/- 8.2 Sv and for heat 2.10 +/- 0.57.The correlation is 0.87.10−5 10−4 10−3 10−2 10−1 100−50−40−30−20−10010EigenValues (Normalized)Magnitude (dB) Spectral Density Relative to Uniform Fourier DistributionEval_funct_U: 05/10/99 10:54:21/home/awi1/mlosch/FLAME/24N/24N_03−Jan−1911/GRACE10/GSerror/isqrSpd680.datChebyshev (40 dB) smoothing, Nspd=680, Nvac=2, Mplt=35010−5 10−4 10−3 10−2 10−1 1000102030405060708090100110Cumulative Spectral Distribution Versus EigenvaluesDistribution function (percent)EigenValues (Normalized)Figure 10:Smoothed spectral density (upper panel) and cumulative spectral distribu-tion (lower panel) of the Hessian of the 24oN problem.Figure 11:Part of the smoothed spectral density weighted for mass (cyan line) and heattransport (purple line) through Florida Strait and its cumulative distribution.Although the dimension of the problem is approximately 12000, only 680coefficients un(x) are needed for the inversion to calculate the posterior errorsof the mass and heat transport through Florida Strait.ConclusionsThree methods are presented for the calculation of posterior errors of inversesolutions. Of these methods, the full inversion of the Hessian is the most accuratebut it is only applicable to small and moderate size problems. Of the other twomethods the conjugate gradients approach leads to reliable results, but is compu-tationally more expensive than finding the solution to the inverse problem itself.The polynomial approximation on the other hand yields also reliable estimates ofthe errors, but it involves an approximation of the Hessian spectrum by Cheby-shev’s Polynomial. This drawback is at the same time a great advantage over theconjugate gradient method because it decreases computing time by an order ofmagnitude.References• Yaremchuk, A. I. and Schröter, J., 1998: Spectral Analysis of Symmetric Operators: Application to the  Laplace Tidal Model. Journal of Computational Physics,147,1-27.• Sloyan, B. M., 1997: The circulation of the Southern Ocean and the adjacent ocean basins determined by inverse methods. Ph.D. thesis. Institute of Antartic andSouthern Ocean Studies, University of Tasmania.• Yaremchuk, M. et. al., 1999: On the zonal and meridional circulation an ocean transports between Tasmania and Antartica, submitted.• Schröter, J., Sloyan,B. M. and Losch, M., 1999: Impact of the GOCE mission for ocean circulation, Midterm report, Task 4.

On the estimation of posterior error bars

Abstract

Similar works

Full text

Available Versions

Electronic Publication Information Center