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

KT K may have similar structure in reconstruction

◮ For an n × n × n image, K has n3 columns! ◮ Noise, e, is white and unknown.

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 3/35

SVD Analysis

Assume we can compute

K = [U1U2] Σ1 Σ2 V T

1

V T

2

• ,

where Σ1 is k × k and corresponds to components such that

uT

i g ≈ uT i ˆ

g.

Noise contaminated exact solution:

f = V1Σ−1

1 (UT 1 g) + dominant

• V2Σ−1

2

(UT

2 g) ≈U T

2 e

Truncated SVD solution: freg = V1Σ−1

1 (UT 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(n3 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−1Kf = M−1g.

If we could compute the SVD of K, the ideal preconditioner would be

M = [U1U2] Σ1 I V T

1

V T

2

• .

Then, semi-convergence in 1 iteration since:

◮ The singular values of M−1K are 1 and Σ2. ◮ M−1g 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

• ur 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      b11C b12C · · · b1nC b21C b22C · · · b2nC . . . . . . . . . . . . bm1C bm2C · · · bmnC      . A ⊗ B ⊗ C is the double-block matrix:      a11B ⊗ C a12B ⊗ C · · · a1nB ⊗ C a21B ⊗ C a22B ⊗ C · · · a2nB ⊗ C . . . . . . . . . . . . am1B ⊗ C am2B ⊗ C · · · amnB ⊗ C      .

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 8/35

SVDs of Kronecker Products

If A = UaΣaV T

a , B = UbΣbV T b , C = UcΣcV T c , then

A ⊗ B ⊗ C = (Ua ⊗ Ub ⊗ Uc)(Σ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 Ai, Bi, Ci with the appropriate structure such that

K −

s

• i=1

Ai ⊗ Bi ⊗ CiF

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 6 1 2 3 7 6 1 2 3 7 6 1 2 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 x h r a b l x h r a l x p q a b x h t p q r a b l x h t p r a l x p q a b t p q r a b t p r a                , P =   h b q x a p l r t  

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 12/35

Triply Toeplitz

Same structure as doubly Toeplitz, except now each Ti 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 ai, bi, ci be the central columns of banded Toeplitz matrices Ai, Bi, Ci, respectively. Then

K −

s

• i=1

Ai ⊗ Bi ⊗ CiF = ¯ P −

s

• i=1

¯ ci ◦ ¯ bi ◦ ¯ aiF,

where the bar notation implies a (diagonal) weighting on the faces of P and the vectors ai, bi, ci, 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

K −

s

• i=1

Ai ⊗ Bi ⊗ CiF = ¯ P −

s

• i=1

¯ ci ◦ ¯ bi ◦ ¯ aiF,

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]

¯ P =

r1

• i=1

r2

• j=1

r3

• k=1

δijkui ◦ vj ◦ wk ≈

s

• m=1

˜ δimuim ◦ vim ◦ wim.

◮ Other, possibilities with little degradation? Work with

Perrone, Martin

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 19/35

Approximation

K −

s

• i=1

Ai ⊗ Bi ⊗ CiF = ¯ P −

s

• i=1

¯ ci ◦ ¯ bi ◦ ¯ aiF

◮ Optimal tensor solution can be calculated if s = 1, used

to construct Ai, Bi, Ci.

◮ 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(n9) flops, storage)

◮ Compute approximation/preconditioner: <= O(n4)

flops

◮ Store a few columns of SVD for 3, n × n matrices.

(minimal)

◮ Preconditioner application O(n3) for small k ◮ Matvec O(n3 log(n))

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 22/35

Example 1

Blurring effects caused by optical limits in 3D microscopy,

• n a stack of 20, 128 × 103 slices of a dendrite. Yields triply

Toeplitz blurring operator. Both examples:

◮ Matlab 7 ◮ RestoreTools [Nagy, Palmer, Perrone, ‘04]

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 23/35

True Image Slices (1,7,13,19)

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 24/35

Blurred Image Slices

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 25/35

Largest 50 |˜

δi|

5 10 15 20 25 30 35 40 45 50 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Magnitude of δ values computed by HOSVD Index of sorted δ values computed by HOSVD

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 26/35

Relative Errors

5 10 15 20 25 30 35 40 45 50 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 Iteration Relative error MRNSD PMRNSD, THOSVD(1) PMRNSD, THOSVD(3)

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 27/35

Reconstructed Images

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 28/35

Example 2

Blurring effects caused by partial volume averaging in spiral CT

κ(x, y, z) = (κ1(x, y, z) + κ2(x, y, z) + κ3(x, y, z))/3 ,

where

κi(x, y, z) = 1

• (2π)3σ3

i

e−(x2+y2+z2)/2σ2

i ,

and σi were chosen randomly with 1 ≤ σi ≤ 2.

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 29/35

True Image Slices 1,8,15,22

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 30/35

Blurred, Noisy Slices 1,8,15,22

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 31/35

Largest |˜

δi|

2 4 6 8 10 12 14 16 18 20 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

Magnitude of δ values computed by HOSVD Index of sorted δ values computed by HOSVD

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 32/35

Convergence History

5 10 15 20 25 30 35 40 45 50 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 Iteration Relative error CGLS PCGLS, THOSVD(1) PCGLS, THOSVD(4)

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 33/35

Restored Images

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 34/35

Conclusions and Future Work

◮ Derived Kronecker product approximations to

• perators in image processing applications, used them

to compute preconditioners making 3D deblurring computationally tractable.

◮ Optimal approximation if one term is used. ◮ Raises interesting questions about tensor

decompositions, which are “best”.

◮ Extension to dense structured matrices (image

processing, PDE’s).

◮ Kronecker approximations do provide useful basis for

edge preserving regularization (work with Per Christian Hansen).

Kronecker Products, Tensor Decompositionsand 3D Imaging Applications – p. 35/35