worklog.org

MAS Worklog

2018-05-05 - CSBS-SSE progress

• rewrote some of the SSE functions to use einsum and more idiomatic python
• decided that sse init() should store the psf ffts for use in later computation
• created a Github repo

SSE Cost for CSBS algorithm

Cost is $\trace (Σ_e) = \trace ((A^HA + λ L^H L)^{-1})$.

We approximate forward model as 2D circular convolutions so that matrices can be diagonalized and we end up with faster computation.

Assume 2 sources.

$$ \begin{aligned} \trace (Σ_e) &= \trace \left ( \begin{bmatrix}F_{2D}^{-1} & 0 \ 0 & F_{2D}^{-1}\end{bmatrix} \begin{bmatrix} Γ_{11} - λ Λ & Γ_{12} \ Γ_{21} & Γ_{22} - λ\Lambda \end{bmatrix}^{-1} \begin{bmatrix}F_{2D} & 0 \ 0 & F_{2D}\end{bmatrix} \right) \ &= \trace \left ( \begin{bmatrix}F_{2D} & 0 \ 0 & F_{2D}\end{bmatrix} \begin{bmatrix}F_{2D}^{-1} & 0 \ 0 & F_{2D}^{-1}\end{bmatrix} \begin{bmatrix} Γ_{11} - λ Λ & Γ_{12} \ Γ_{21} & Γ_{22} - λ\Lambda \end{bmatrix}^{-1} \right) \ &= \trace \left ( \begin{bmatrix} Γ_{11} - λ Λ & Γ_{12} \ Γ_{21} & Γ_{22} - λ\Lambda \end{bmatrix}^{-1} \right) \end{aligned} $$

Outline of CSBS Algorithm

• Initialization

  Intialization code before CSBS begins

• Iterate

  Repeat N times

  • for each row of $A$
    • temporarily remove row from $A$ and compute cost of this reduced $A$
  • permanently remove whichever row of $A$ incurred the lowest cost

Note: 1 row of A is a row vector of convolution matrices, length $S$

CSBS Optimizations for SSE

We make a number of optimizations when calculating the SSE cost to make it computationally feasible.

• store “compressed” form of diagonal matrices, PSFs

  $Γ_{ij}$ and $Λ$ are diagonal matrices of dimension $(m ⋅ n) × (m ⋅ n)$. we store the diagonal elements only in a matrix of size $m × n$

  We store the PSFs (and DFTs) in their 2D form to avoid flattening/reshaping. Consequency, we have special multiplication and inversion functions which operate on these compressed matrices directly.

• update $Γ = \begin{bmatrix} Γ_{11} & Γ_{12} \ Γ_{21} & Γ_{22} \end{bmatrix}$ instead of recomputing it

  $Γ$ changes very little between CSBS iterations. We subtract off the contribution from the removed row of $A$ at the end of each CSBS iteration

• dont store duplicate rows in $A$

  $A$ contains many duplicate rows to represent repeated measurements at one measurement plane. We only store one copy of a row in $A$ along with a counter of how many copies of this row are left

• use block matrix inversion formula for inverting $\bar{Γ} = \begin{bmatrix} Γ_{11} - λ Λ & Γ_{12} \ Γ_{21} Λ & Γ_{22} - λ \end{bmatrix}$

  The elements of $\bar{Γ}$ are diagonal matrices. We use the block matrix inversion formula to invert it efficiently.

Outline of CSBS Algorithm for SSE

• Initialization
  • precompute PSF ffts
  • calculate full $Γ$ matrix and $Λ$
• Iterate
  • for each row of $A$
    • temporarily remove row from $A$ and compute cost: $\trace (Σ_e)$
  • permanently remove whichever row of $A$ incurred the lowest cost
  • update $Γ$ by subtracting contribution from removed row

2018-05-16 - CSBS-SSE completion

• fixed block_inv to work with even sized inputs
• added iteration_end function to cost_module
• derived $Dx^T ⋅ Dx + Dy^T ⋅ Dy$

Finding $Λ$

$D_x^T D_x + D_y^T D_y$ is a block circulant matrix with circulant blocks, so it can be diagonalized by the 2D DFT matrix, $F_{2D}$

$$Λ = F_{2D} (D_x^T D_x + D_y^T D_y) F_{2D}^{-1}$$

Then the compressed form of $Λ$ is

$$\text{vect}^{-1}(\text{diag}(Λ)) = \text{2D DFT of vect}^{-1}(\text{1st row of } D_x^T D_x + D_y^T D_y)$$

$D_x$ and $D_y$ are discrete derivative operators operating on a flattened image $x$ of size $m × n$

$$D_x x = \text{vect}\left(\begin{bmatrix} \horzbar & d \circledast x_{r1} & \horzbar & \vdots & \\ \horzbar & d \circledast x_{rm} & \horzbar \\ \end{bmatrix}\right) = \text{vect}\left(\begin{bmatrix} \horzbar & D_r x_{r1} & \horzbar \\ & \vdots & \\ \horzbar & D_r x_{rm} & \horzbar \\ \end{bmatrix}\right) $$

where $D_r$ is the 1D circulant matrix of $[-1, 1]$ of size $n× n$. $\text{vect}$ is an operator which concatenates the rows of a matrix into a single vector.

Similarly for $D_y$

$$D_y x = \text{vect}\left(\begin{bmatrix} \vertbar & & \vertbar d \circledast x_{c1} & \hdots & d \circledast x_{cn} \\ \vertbar & & \vertbar \\ \end{bmatrix}\right) = \text{vect}\left(\begin{bmatrix} \vertbar & & \vertbar \\ D_c x_{c1} & \hdots & D_c x_{cn} \\ \vertbar & & \vertbar \\ \end{bmatrix}\right) $$

where $D_c$ is the 1D circulant matrix of $[-1, 1]$ of size $m× m$

$$D_x = I_{m × m} ⊗ D_r$$ $$D_y = D_c ⊗ I_{n × n}$$

$$D_x^T D_x + D_y^T D_y = (I_{m × m} ⊗ D_r)^T (I_{m × m} ⊗ D_r) + (D_c ⊗ I_{n × n})^T (D_c ⊗ I_{n × n})$$

Using the property $(A ⊗ B)(C ⊗ D) = AC ⊗ BD$, we get

$$D_x^T D_x + D_y^T D_y = I_{m × m} ⊗ D_r^TD_r + D_c^TD_c ⊗ I_{n × n}$$

Substituting back in,

$$ \begin{aligned} & \text{2D DFT of vect}^{-1}(\text{1st row of } D_x^T D_x + D_y^T D_y) &= \text{2D DFT of vect}^{-1}(\text{1st row of } I_{m × m} ⊗ D_r^TD_r + D_c^TD_c ⊗ I_{n × n}) \\ &= \text{2D DFT of } \left[ \text{vect}^{-1}(\text{1st row of } I_{m × m} ⊗ D_r^TD_r) + \text{vect}^{-1}(\text{1st row of } D_c^TD_c ⊗ I_{n × n}) \right] \end{aligned}$$

---

Let $A$ and $B$ be arbitrary matrices

$$I ⊗ A = \begin{bmatrix} A & & \ & \ddots & \ & & A \end{bmatrix}$$

$$\text{vect}^{-1}(\text{1st row of } I ⊗ A) = \begin{bmatrix} a_{11} & \hdots & a_{1n} \ 0 & \hdots & 0 \ \vdots & & \vdots \ 0 & \hdots & 0 \end{bmatrix}$$

$$B ⊗ I = \begin{bmatrix} b_{11} \ident & \hdots & b_{1n} \ident & \vdots & \\ b_{m1} \ident & \hdots & b_{mn} \ident \end{bmatrix}$$

$$\text{vect}^{-1}(\text{1st row of } B ⊗ I) = \begin{bmatrix} b_{11} & 0 & \hdots & 0 \ \vdots & \vdots & & \vdots \ b_{1n} & 0 & \hdots & 0 \end{bmatrix}$$

---

If we let $A = D_x^TD_x$ and $B = D_y^TD_y$, then

$$\text{2D DFT of } \left[ \text{vect}^{-1}(\text{1st row of } I_{m × m} ⊗ D_r^TD_r) + \text{vect}^{-1}(\text{1st row of } D_c^TD_c ⊗ I_{n × n}) \right] = \begin{bmatrix} a_{11} + b_{11} & a_{12} & \hdots & a_{1n} b_{12} & 0 & \hdots & 0 \\ \vdots & \vdots & & \vdots \\ b_{1n} & 0 & \hdots & 0 \end{bmatrix}$$

2018-05-20 - Framework

Today we began work on a mathematical framework to formalize the constraints and goals of the plane selection/exposure time problem.

Parameters

Parameters we control in the problem.

• exposure time $τ_k$
• measurement plane locations $d_k$
• measurement plane transition time $Δ_j$

Goals

Problem optimization goals.

• high SNR (maximize $τ_k$)
• Minimize temporal artifacts (minimize $τ_k$, minimize $Δ_j$)
• Capture measurement diversity (maximize order of $d_k$)

3 types of noise

             --------------       -----------------------------
source---+---| microphone |-------| system processing ----+---|-------
         |   --------------       ------------------------|----
        n_2                                              n_3

$y = n_1(Ax) + A n_2 + n_3$

• $n(Ax)$ - shot noise. large input signal increases self interference
• $n_2$ - dark noise (environmental noise). e.g. computer fan
• $n_3$ - read noise (system noise). e.g. ADC noise, self interference

2018-05-22 - Time considerations of the PSSI

We are trying to image a dynamically changing object. Hence, we cannot keep exposure times very long. We also need to consider the transition time of the detector between measurement planes. Here, we formulate these, and set a condition to satisfy:

Parameters

• number of measurement planes : $K$
• exposure time at each measurement plane : $t_{exp}$
• detector transition time from $i^{th}$ to $(i+1)^{th}$ measurement plane: $t_{tr}^{i}$
• the time for which the dynamic object can be considered still: $t_{obj}$

Requirement

The total time to complete taking measurements should not exceed $t_{obj}$:

• $K t_{exp} + ∑_{i=1}^{K-1} t_{tr}^{i} \leq t_{obj}$

2018-05-22 - Plotting CSBS Results

Wrote some code to visualize output from CSBS, shown below.

I noticed that with a poor choice of lambda, CSBS sometimes completely kills off other focused measurement planes.