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
Recommend
More recommend
