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

Kronecker Products, Tensor Decompositions and 3D Imaging Applications Misha E. Kilmer Tufts University James Nagy Emory University Lisa Perrone Hawaii Pacific University Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 1/35

Outline ◮ Preconditioning for discrete ill-posed problems ◮ The matrix approximation problem ◮ The role of tensors ◮ Theoretical Results ◮ Numerical results ◮ Conclusions and future work Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 2/35

Problem Description Model: Kf = ˆ g + e = g ◮ K is ill-conditioned, no gap in SV spectrum ◮ K is triply Toeplitz or triply T+H 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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 3/35

SVD Analysis Assume 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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 4/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 5/35

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 so no noise introduced. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 6/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 7/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 8/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 9/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 10/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 11/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 12/35

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 3rd order tensor P . Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 13/35

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. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 14/35

Matrix Approximation Problem Revisited Note that if A, B, C are banded Toeplitz (or T+H), they are uniquely specified by their respective central columns, vectors a, b, c . Also, A ⊗ B ⊗ C is triply Toeplitz (T+H) with special banding structure. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 15/35

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 (diagonal) weighting on the faces of P and the vectors a i , b i , c i , and ◦ implies 3way outer product. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 16/35

Theorem 2 Similar result when K is triply Toeplitz + Hankel (banded), except the weighting matrix is a very special upper triangular matrix. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 17/35

Tensor approximation s s � � c i ◦ ¯ A i ⊗ B i ⊗ C i � F = � ¯ � K − P − ¯ b i ◦ ¯ a i � F , i =1 i =1 Computing optimal approximation to K requires computing the optimal rank-s approximation to the tensor ¯ P . ◮ For s = 1 , unique solution. ◮ For s > 1 , is there a solution? (uniqueness requires mild constaints). Can we compute it? ◮ A good choice for small s may not be known a priori. ◮ There are no “orthogonality” constraints on the decomposition, and none needed for our application. Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 18/35

Tensor Decompositions Choices: ◮ PARAFAC model (N-way Toolbox 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 . i =1 j =1 m =1 k =1 ◮ Other, possibilities with little degradation? Work with Perrone, Martin Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 19/35

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 can be 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. ◮ The HOSVD for the tensor should not be confused with the SVD of K . Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 20/35

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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 21/35

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 4 ) 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 Products, Tensor Decompositionsand 3D Imaging Applications – p. 22/35