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)

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 2/33

SVD Analysis

Assume, temporarily, 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 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(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 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−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 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

• 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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 6/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 7/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 8/33

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 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 6 1 2 3 7 6 1 2 3 7 6 1 2 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 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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 11/33

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 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 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 weighting on the faces of P and the vectors ai, bi, ci, 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 Ai ⊗ Bi ⊗ Ci requires we be able to compute the

• ptimal 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]

¯ P =

r1

• i=1

r2

• j=1

r3

• k=1

δijkui ◦ vj ◦ wk ≈

s

• m=1

˜ δimuim ◦ vim ◦ wim.

◮ We are working on other possibilities...

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 17/33

Approximation

K −

s

• i=1

Ai ⊗ Bi ⊗ CiF = ¯ P −

s

• i=1

¯ ci ◦ ¯ bi ◦ ¯ aiF

◮ Optimal tensor solution 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.

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

◮ Compute approximation/preconditioner: O(n3) 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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 20/33

Example 1

Blurring effects, caused by optical limits in 3D microscopy,

• n 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|

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 24/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 25/33

Reconstructed Images

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 26/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 27/33

True Image Slices 1,8,15,22

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 28/33

Blurred, Noisy Slices 1,8,15,22

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 29/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 30/33

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 Product Approximation for Preconditioningin 3D Imaging Applications – p. 31/33

Restored Images

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 32/33

Conclusions and Future Work

◮ Derived Kronecker product approximations to

• perators in image processing applications, used them

to compute preconditioners making 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 may provide useful basis for

edge preserving regularization (work with Per Christian Hansen).

Kronecker Product Approximation for Preconditioningin 3D Imaging Applications – p. 33/33