Randomized Matrix-Free Trace and Log-Determinant Estimators

31 Pages·2017·3.28 MB·English

Randomized Matrix-Free
Trace and Log-Determinant Estimators
Ilse Ipsen
Joint with: Alen Alexanderian & Arvind Saibaba
NorthCarolinaStateUniversity
Raleigh,NC,USA
ResearchsupportedbyDARPAXData
Our Contribution to this Minisymposium
Inverse problems: Bayesian OED
Big data: Randomized estimators
Given: Hermitian positive semi-definite matrix A
Want: trace(A) and logdet(I+A)
Our estimators
• Matrix free, simple implementation
• Much higher accuracy than Monte Carlo
• Informative probabilistic bounds, even for small dimensions
Inverse Problem: Diffusive Contaminant Transport
Forwardproblem: Time-dependentadvection-diffusionequation
Inverseproblem: Re-constructuncertaininitialconcentration
frommeasurementsof35sensors,at3timepoints
Prior-preconditionedFisherinformation H≡C1/2F∗Γ−1 FC1/2
prior noise prior
SensitivityofOED: trace(H)
BayesianD-optimaldesigncriterion: logdet(I+H)
Our estimators
• Fast and accurate for Bayesian inverse problem
• Matrix free, simple implementation
• Much higher accuracy than Monte Carlo
• Informative probabilistic bounds, even for small dimensions
Our Contribution to this Minisymposium
Inverse problems: Bayesian OED
Big data: Randomized estimators
Given: Hermitian positive semi-definite matrix A
Want: trace(A) and logdet(I+A)
Our Contribution to this Minisymposium
Inverse problems: Bayesian OED
Big data: Randomized estimators
Given: Hermitian positive semi-definite matrix A
Want: trace(A) and logdet(I+A)
Our estimators
• Fast and accurate for Bayesian inverse problem
• Matrix free, simple implementation
• Much higher accuracy than Monte Carlo
• Informative probabilistic bounds, even for small dimensions
Existing Estimators
Small Explicit Matrices: Direct Methods
Trace = sum of diagonal elements
(cid:88)
trace(A) = a
jj
j
Logdet
(1) Convert to trace
logdet(I+A) = trace(log(I+A))
(2) Cholesky factorization I+A = LL∗
(cid:89)
logdet(I+A) = log|det(L)|2 = 2log l
jj
j
Large Sparse or Implicit Matrices
Trace: Monte Carlo methods
N
(cid:88)
trace(A) ≈ 1 z∗Az
N j j
j=1
N independent random vectors z
j
Rademacher, standard Gaussian ⇒ unbiased estimator
Hutchinson1989, Avron&Toledo2011
Roosta-Khorasani&Ascher2015, Lin2016
Logdet: Expansion of log
 
(cid:88)
logdet(I+A) ≈ trace pj(A)
j
Barry&Pace1999, Pace&LeSage2004, Zhangetal. 2008
Chenetal. 2011, Anitescuetal. 2012,Boutsidisetal. 2015
Hanetal. 2015
Our Estimators
Randomized Subspace Iteration
Given: n×n matrix A with k dominant eigenvalues
Compute low rank approximation T:
(1) Pick random starting guess Ω with (cid:96) ≥ k columns
(2) Subspace iteration Y = AqΩ
(3) Orthonormalize: Thin QR Y = QR
(4) Compute (cid:96)×(cid:96) matrix T = Q∗AQ
Estimators:
trace(T) ≈ trace(A) logdet(I+T) ≈ logdet(I+A)