Misha E. Kilmer Tufts University James Nagy Emory University Lisa - PowerPoint PPT Presentation

Kronecker Product Approximation for Preconditioning in 3D Imaging Applications Misha E. Kilmer Tufts University James Nagy Emory University Lisa Perrone Tufts University Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 1/33

Problem Description Model: Kf = ˆ g + e = g ◮ K is ill-conditioned, no gap in SV spectrum ◮ K is structured: * K is triply Toeplitz or T+H (restoration) K T K may have similar structure in reconstruction ◮ For an n × n × n image, K has n 3 columns! ◮ Noise, e , is white and unknown. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 2/33

SVD Analysis Assume, temporarily, we can compute � Σ 1 � � V T � 0 1 K = [ U 1 U 2 ] , V T 0 Σ 2 2 where Σ 1 is k × k and corresponds to components such that u T i g ≈ u T i ˆ g. Noise contaminated exact solution: dominant � �� � f = V 1 Σ − 1 1 ( U T V 2 Σ − 1 ( U T 1 g ) + 2 g ) 2 � �� � ≈ U T 2 e Truncated SVD solution: f reg = V 1 Σ − 1 1 ( U T 1 g ) Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 3/33

Iterative Regularization We cannot compute the SVD of K ! But K is structured, so we can compute matrix-vector products using FFT’s quickly O ( n 3 log n ) flops if image is n × n × n . This means we should use an iterative regularization scheme (eg. CGLS, MRNSD). We stop iterating before the solution converges to the exact solution of the system. Cost is O ( Nj ) where a matvec costs O ( N ) and j is the number of iterations until semi-convergence. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 4/33

Ideal Preconditioning Major difficulty is that j can be large! Consider the left preconditioned system M − 1 Kf = M − 1 g. If we could compute the SVD of K , the ideal preconditioner would be � Σ 1 � � V T � 0 1 M = [ U 1 U 2 ] . V T 0 I 2 Then, semi-convergence in 1 iteration since: ◮ The singular values of M − 1 K are 1 and Σ 2 . ◮ M − 1 g looks similar to the TSVD solution and therefore no noise is introduced. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 5/33

Goal Since we cannot compute the SVD of K , the main goal of our work is to approximate K by a matrix for which we can compute (efficiently) and store (implicitly) the SVD, and use the SVD of the approximation to construct M . Key : Using the structure of K , the matrix approximation problem can be reduced to a computationally tractable problem involving 3-way arrays (tensors). Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 6/33

Kronecker Products B ⊗ C is the block matrix   b 11 C b 12 C · · · b 1 n C b 21 C b 22 C · · · b 2 n C      . . . . .   . . . . . . . .  b m 1 C b m 2 C · · · b mn C A ⊗ B ⊗ C is the double-block matrix:   a 11 B ⊗ C a 12 B ⊗ C · · · a 1 n B ⊗ C a 21 B ⊗ C a 22 B ⊗ C · · · a 2 n B ⊗ C      . . . . .   . . . . . . . .  a m 1 B ⊗ C a m 2 B ⊗ C · · · a mn B ⊗ C Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 7/33

SVDs of Kronecker Products If A = U a Σ a V T a , B = U b Σ b V T b , C = U c Σ c V T c , then A ⊗ B ⊗ C = ( U a ⊗ U b ⊗ U c )(Σ a ⊗ Σ b ⊗ Σ c )( V T a ⊗ V T b ⊗ V T c ) , which is an SVD (under appropriate ordering). Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 8/33

Matrix Approximation Problem Find A i , B i , C i with the appropriate structure such that s � � K − A i ⊗ B i ⊗ C i � F i =1 is minimized. (reason for s > 1 will be discussed later) Related 2D work: ◮ Kamm & Nagy, ’00 ◮ Nagy, & Ng, Perrone, ’04 ◮ Perrone, ’05 ◮ K. & Nagy, ’06 Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 9/33

Banded Toeplitz Complete information about this banded Toeplitz matrix is captured in a central vector of length n :   1 2 3 0 0   6 1 2 3 0     7 6 1 2 3 .     0 7 6 1 2   0 0 7 6 1 Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 10/33

Doubly Toeplitz This doubly Toeplitz matrix can be represented by its central column, which we reshape:   a b 0 x h 0 0 0 0 r a b l x h 0 0 0     0 r a 0 l x 0 0 0         p q 0 a b 0 x h 0 h b q     t p q r a b l x h , P = x a p       0 t p 0 r a 0 l x l r t     0 0 0 p q 0 a b 0     0 0 0 t p q r a b   0 0 0 0 t p 0 r a Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 11/33

Triply Toeplitz Same structure as doubly Toeplitz, except now each T i is a doubly Toeplitz matrix. Similarly, if banded structure, the matrix is completely represented by its central column , which reshape into a 3-way array (tensor) P . Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 12/33

Structure If we assume that the blurring in the 3D image is spatially invariant and that the image satisfies 0 boundary conditions, K will be triply Toeplitz with this banded structure. If we assume reflexive boundary conditions, K will be triply Toeplitz+Hankel with special banding: For ease of discussion, we omit this case. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 13/33

Matrix Approximation Problem Revisited Note that if A, B, C are banded Toeplitz, they are uniquely specified by their respective central columns, vectors a, b, c . Also, A ⊗ B ⊗ C is triply Toeplitz (banded). Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 14/33

Theorem 1 Let K be triply Toeplitz (banded) and P be the 3D tensor defining the central column of K . Let a i , b i , c i be the central columns of banded Toeplitz matrices A i , B i , C i , respectively. Then s s � � A i ⊗ B i ⊗ C i � F = � ¯ c i ◦ ¯ � K − P − ¯ b i ◦ ¯ a i � F , i =1 i =1 where the bar notation implies a weighting on the faces of P and the vectors a i , b i , c i , and ◦ implies 3way outer product. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 15/33

Tensor approximation Being able to compute the optimal approximation to K by � s i =1 A i ⊗ B i ⊗ C i requires we be able to compute the optimal rank-s approximation to the tensor ¯ P . ◮ For s = 1 , unique solution. ◮ For s > 1 , is there a unique solution? Can we compute it? ◮ There are no “orthogonality” constraints on the decomposition. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 16/33

Tensor Decompositions Choices: ◮ PARAFAC model (package by R. Bro) ◮ HOSVD [ de Lathauwer, de Moor and Vandewalle, ‘00] r 1 r 2 r 3 s � � � � ˜ ¯ P = δ ijk u i ◦ v j ◦ w k ≈ δ i m u i m ◦ v i m ◦ w i m . m =1 i =1 j =1 k =1 ◮ We are working on other possibilities... Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 17/33

Approximation s s � � c i ◦ ¯ A i ⊗ B i ⊗ C i � F = � ¯ � K − P − ¯ b i ◦ ¯ a i � F i =1 i =1 ◮ Optimal tensor solution calculated if s = 1 , used to construct A i , B i , C i . ◮ Settle for HOSVD approximation for s > 1 . ◮ To determine s , observe | ˜ δ i | . Gaps/decay for small values of s imply K is well approximated by small number of terms in the Kronecker sum. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 18/33

Preconditioner ◮ Case s = 1 , compute SVD from Kronecker product approximation directly and use it to construct M . ◮ Case s > 1 , one more level of approximation required (skip details), but M still ends up specified in Kronecker form. Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 19/33

Compute/Storage Summary For an n × n × n image (an SVD of K would cost O ( n 9 ) flops, storage) ◮ Compute approximation/preconditioner: O ( n 3 ) flops ◮ Store a few columns of SVD for 3, n × n matrices. (minimal) ◮ Preconditioner application O ( n 3 ) for small k ◮ Matvec O ( n 3 log( n )) Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 20/33

Example 1 Blurring effects, caused by optical limits in 3D microscopy, on a stack of 20 , 128 × 103 slices of a dendrite. Both examples: ◮ Matlab 7 ◮ RestoreTools [Nagy, Palmer, Perrone, ‘04] Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 21/33

True Image Slices (1,7,13,19) Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 22/33

Blurred Image Slices Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 23/33

Largest 50 | ˜ δ i | 1 10 0 10 Magnitude of δ values computed by HOSVD −1 10 −2 10 −3 10 −4 10 −5 10 −6 10 0 5 10 15 20 25 30 35 40 45 50 Index of sorted δ values computed by HOSVD Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 24/33