Note

This page was generated from docs/notebooks/inversion/rtm2svd.ipynb.

Compute SVD from Raytransfer Matrix#

When it comes to tomographic inversion algorithms, SVD calculation is often needed and cherab.phix.inversion subpackage also requires. An RTM is a so large matrix that we have to prepare a lot of computer resources with large memory size.

SVD is performed by the following equation:

\[U\Sigma V^\mathsf{T} = AL^{-1},\]

where \(A\) is an RTM, \(L\) is the regularization operator (e.g. \(L=\text{Laplacian}\) if using Phillips regularization).

We proceed with the following process to compute SVD:

  1. Compute \(L^{-1}\)

  2. Compute the product of matrices \(AL^{-1}\)

  3. Compute SVD with \(AL^{-1}\)

To reduce memory usage, each generated arrays are stored on disk and deleted from RAM. Additionally, we compute the folloing matrix for the future inversion work.

  1. Compute \(L^{-1}V\)

[1]:

import os
import platform
from pathlib import Path
from textwrap import dedent

import numpy as np
import psutil
from bokeh.plotting import output_notebook

# To show resource usage during calculation
from dask.diagnostics import ResourceProfiler
from numpy.lib.format import open_memmap
from raysect.optical import World
from scipy.sparse.linalg import inv

from cherab.phix.tools import laplacian_matrix
from cherab.phix.tools.raytransfer import import_phix_rtc

output_notebook()

# directory path where the RTM is stored.
RTM_DIR = Path().cwd().parent.parent.parent / "output" / "RTM" / "2022_12_13_00_49_29"

# Load PHiX Raytransfer object
world = World()
rtc = import_phix_rtc(world)

# show system info
core = os.cpu_count()
gb = float(1024**3)
mem_total = int(psutil.virtual_memory()[0] / gb)
storage_total = int(psutil.disk_usage("/")[0] / gb)
print(
    dedent(
        f"""
        ---------- System Info ----------
        System   : {platform.system()}, {platform.machine()}
        Kernel   : {platform.release()}
        Compiler : {platform.python_compiler()}
        CPU      : {core} (Core)
        Memory   : {mem_total} GiB
        Disk     : {storage_total} GiB
        ---------------------------------
        """
    )[1:-1]
)
Loading BokehJS ... 
---------- System Info ----------
System   : Linux, x86_64
Kernel   : 6.2.0-34-generic
Compiler : GCC 11.3.0
CPU      : 128 (Core)
Memory   : 125 GiB
Disk     : 915 GiB
---------------------------------

Compute \(L^{-1}\)#

[2]:

# define resource profiler to monitor usage of memory and CPU
resource_profiler = ResourceProfiler()

with resource_profiler:
    # Compute Laplacian sparse matrix
    laplacian = laplacian_matrix(rtc.voxel_map, dir=8)

    # create memory-map to store array
    L_inv = open_memmap(
        RTM_DIR / "L_inv.npy", dtype=np.float64, mode="w+", shape=(rtc.bins, rtc.bins)
    )

    # Compute L^-1
    L_inv[:] = inv(laplacian.tocsc()).toarray()

resource_profiler.visualize()

[2]:
figure(
id = 'p1001', …)

Compute \(AL^{-1}\)#

Load RTM with the reshape to 2D array from 3D

[3]:

rtm = open_memmap(RTM_DIR / "rtm.npy").reshape((-1, rtc.bins))

[4]:

# create memory-map to store array
AL_inv = open_memmap(RTM_DIR / "AL_inv.npy", dtype=np.float64, mode="w+", shape=rtm.shape)

# compute AL^-1
with resource_profiler:
    AL_inv[:] = rtm.dot(L_inv)

resource_profiler.visualize()

[4]:
figure(
id = 'p1057', …)

Compute SVD: \(U\Sigma V^\mathsf{T} = AL^{-1}\)#

Define \(U, \Sigma, V^\mathsf{T}\) memory-map array firstly.

[5]:

# define array shape number
m, n = AL_inv.shape
k = min(n, m)

u_row = k if m < n else m
vh_col = n if m < n else k

# create memory-map to store
u = open_memmap(RTM_DIR / "u.npy", dtype=np.float64, mode="w+", shape=(u_row, k))
s = open_memmap(RTM_DIR / "s.npy", dtype=np.float64, mode="w+", shape=(k,))
vh = open_memmap(RTM_DIR / "vh.npy", dtype=np.float64, mode="w+", shape=(k, vh_col))

compute SVD

[6]:

# compute SVD
with resource_profiler:
    u[:], s[:], vh[:] = np.linalg.svd(AL_inv, full_matrices=False)

resource_profiler.visualize()

[6]:
figure(
id = 'p1113', …)

Compute \(L^{-1}V\)#

[7]:

# create memory-map to store
L_inv_V = open_memmap(RTM_DIR / "L_inv_V.npy", dtype=np.float64, mode="w+", shape=(n, k))

# compute
with resource_profiler:
    L_inv_V[:] = L_inv.dot(vh.T)

resource_profiler.visualize()

[7]:
figure(
id = 'p1169', …)

Note

demos/rtm2svd.py is the script file having the same workflow of SVD calculation as the above codes.