SVDによる圧縮の利点と限界点について

挨拶

フリュー(株)で、プリントシール機のソフト開発をしている金(キム)です。 よろしくお願いいたします！

序論

目的

本記事の目的は、特異値分解（Singular Value Decomposition, SVD）を用いて画像圧縮を行い、その有用性と実用性を検討することです。SVDによる低ランク近似を活用し、圧縮率と画質のバランスを評価しました。また、JpegやWebPと比較しながら、SVDを画像圧縮に適用する際の利点や課題を明らかにします。

前提知識

特異値分解の計算の流れを説明するため以下の知識が必要となります。
1．線形代数の理論（行列演算、直行行列＆対角行列、固有値＆固有ベクトル、主成分分析 ）
2．整数と実数の情報に必要なbitのサイズ
3．整数型と実数型の表現可能範囲
4．画質や評価方法（PSNR、SSIM、MSE）
5．pythonの簡単な文法

SVDの理論

SVD

特異値分解（Singular Value Decomposition, SVD）は、任意の A を、三つの特定の行列の積に分解する手法です。各行列は、行列によって変換された後も、その方向が変わらない固有ベクトルとその方向の大きさを表す固有値で構成され、いくつかの特徴を持っています。UはAの固有ベクトルで式が誘導されるため行列Aの特徴を表し、Vᵀは行空間として再構成された固有ベクトルで表現されるためAの新しい基底での表現を示します。そして、Σは固有値を平方根して得た特異値で表現されるためAの各方向の大きさを示しています。このような特徴を持つSVDについて、具体的な内容を整理すると、次のようになります。

A：元の行列（m×nのサイズ）
U：m×mの直交行列で、左特異ベクトルを含む
Σ：m×nの対角行列で、特異値が対角要素として配置される
Vᵀ：n×nの直交行列で、右特異ベクトルを含む

UとVの特徴

ここでUとVは直交行列です。直交行列とは、行または列が互いに直交しており、各ベクトルの大きさが1である行列です。直交行列の特徴は、転置行列を掛けると単位行列になることです。数学的には以下のように表されます。

UUᵀ=UᵀU=I 、VVᵀ=VᵀV=I

Uは前述のように固有ベクトルの特性を持つため入力空間を回転または反射させる役割を持っています(行列Aがどんな方向から拡張されたのか、その特徴=方向を表す)。Uは元の行列 Aの左特異ベクトル（left singular vectors）を含み、このベクトルは行列Aの列空間を構成します。ここで列空間は元の行列Aの列要素がどの方向に広がっているか（列空間の基底）を表します。

Vも同じく固有ベクトルでの特性を持つため出力空間を回転または反射させる役割を持っています。Vは元の行列Aの右特異ベクトル（right singular vectors）を含み、このベクトルは行列Aの行空間を構成します。ここで行空間は元の行列Aの行の要素がどの方向に広がっているか（行空間の基底）を表します。
回転と反射は、線形変換において重要な役割を果たします。UとVによって元の空間を単純に回転させたり角度を変えることで、元の行列の構造的特性を表現することができます。ただし、Vは行空間で固有ベクトルが新しく解析されたことに気を付ける必要があります。

Σの特徴

Σは対角行列で、主対角成分に特異値（singular values）が並んでいます。

ここで σ₁,σ₂,σ₃,… は行列 A の特異値です。これらの特異値はすべて0以上の実数であり、通常は降順に並べられます（ σ₁ ≥ σ₂ ≥ σ₃ ≥⋯≥0）。
Σ は拡大縮小の役割を持っていて主対角線の値が各方向のベクトルの大きさを縮小または拡大するのに使用されます。大きな特異値は対応する方向で大きさが大きく拡張されていることを意味し、小さな特異値は縮小されていることを意味します。これにより、元の行列Aがどの軸方向でどれだけ強く影響を与えているのかがわかります。

固有ベクトルが行列の特徴を表す理由

固有ベクトルは、行列Aに対して次の式を満たすベクトルvです。

Av=λv

ここで、λは対応する固有値と呼ばれる数値です。この式は、行列Aによる線形変換が、固有ベクトルvの方向を変えず、λによってその大きさだけを変えることを意味します。 固有ベクトルは行列を表すため必要な主軸です。例えば、以下のように固有ベクトルが２つ存在すると２つの固有ベクトルの線形結合で元の行列を表現することができます。行列を表現する軸になりますので、行列の特徴を表していると考えれらます。

SVDを用いた画像圧縮

近似行列A₁

A を構成する項の中で、σ₁、u₁、vᵀ₁はAの最も重要な部分（最大の寄与）を表します。具体的には、
1．σ₁はAの中で最も大きな特異値（Aの主な力）
2．u₁はAの列空間で最も重要な方向を表す単位ベクトル
3．vᵀ₁はAの行空間で最も重要な方向を表す単位ベクトル
を表しています。
したがって、A₁=σ₁u₁vᵀ₁は元のAの中で最も情報量が多い方向のベクトルを抽出していると考えられます。元の行列の主成分や特徴量の抽出に使われます。

画像圧縮の原理と条件

①ピクセルの単純化

まず、画像の各チャンネルを分けてSVDを適用することでピクセルを単純化することができます。以下の図で画像の一部をSVDしてから特異値を削った結果画像のピクセルが元より単純化されたのが分かります。この単純化された行列を画像形式で保存すると、近辺の要素の値が一致した場合、画像を効率よく保存するため画像の容量が減ります。SVDによる画像圧縮は特に専用の形式が存在しないため、他の画像形式を使うか独自で作るしかないですが、画像を構成するピクセルを単純化することで圧縮が効率よくされるようになります。

②行列の要素を削減

画像を読み込み、行列を扱う際にメモリ上の占有容量を減らしたい場合、SVDが活用されられます。基本的な方法として、①と同じく特異値を選定し小さい行列を使うことで、要素数を減らせばメモリ上に格納される情報の容量を減らせます。最初SVDを用いて分解された各行列の総要素数は実は元の行列の要素数より多くなるため（近似行列要素数：m x k + k + n x k）、以下の条件を満足する特異値の数kを厳選する必要があります。

1．行列Aの要素が4BytesでUΣVᵀの要素が4Bytesの場合：m x n > k(m+n+1) → (m x n) / (m + n + 1) > k
2．行列Aの要素が1Bytes(例：0~255)でUΣVᵀの要素が4Bytesの場合： m x n > 4k(m+n+1) → (m x n) / 4(m + n + 1) > k

ここで注意しなければならない所は、画像によっては2Bytesや4Bytesで表現される（医療画像など）こともありますが、一般的には各要素が大体1Byte（0~255）で表現されることです。SVD適用後は各行列の要素が小数点を含む実数情報（情報損失を防ぐため4Bytes）になるため、元の画像で一つの要素数を表現するために必要な容量を確認し、SVD適用後のすべての要素が必要としている容量を計算する必要があります。

SVDを活用した画像圧縮方法

計算の流れは以下になります。pythonスクリプトを追加しましたので、ご参考ください。
①行列の定義と表示
②転置行列の計算
③AᵀA を計算
④固有値の計算
⑤固有ベクトルの計算
⑥特異値の計算(√λ (λ=固有値))
⑦行列 V の構築
⑧行列 U の計算(A=UΣVᵀ -> U=AVΣ⁻¹またはAAᵀの固有ベクトルで構築)
⑨行列 A の再構築(A=UΣVᵀ)
⑩画像圧縮

スクリプトを確認する

①行列の定義と表示

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# 画像を読み込む（小さなサイズの画像を使用）
image_path = '/content/image.jpg'  # 実際の画像ファイルのパスに置き換えてください
image = Image.open(image_path)
image = np.array(image)

# RGBチャンネルに分割
R, G, B = image[:,:,0], image[:,:,1], image[:,:,2]

# 画像をNumPy配列に変換
A_r = np.array(R, dtype=uint8)
A_g = np.array(G, dtype=uint8)
A_b = np.array(B, dtype=uint8)

# 画像を表示
fig, axs = plt.subplots(1, 4, figsize=(10, 5))

axs[0].imshow(image)
axs[0].set_title('Original image')
axs[0].axis('off')

axs[1].imshow(A_r)
axs[1].set_title('R channel')
axs[1].axis('off')

axs[2].imshow(A_b)
axs[2].set_title('G channel')
axs[2].axis('off')

axs[3].imshow(A_g)
axs[3].set_title('B channel')
axs[3].axis('off')
plt.show()

②転置行列の計算

# 行列の転置を計算
A_rT = A_r.T
A_gT = A_g.T
A_bT = A_b.T

fig, axs = plt.subplots(1, 3, figsize=(10, 5))

# 転置行列を表示
axs[0].imshow(A_rT)
axs[0].set_title('Transposed Image(R)')
axs[0].axis('off')

axs[1].imshow(A_gT)
axs[1].set_title('Transposed Image(G)')
axs[1].axis('off')

axs[2].imshow(A_bT)
axs[2].set_title('Transposed Image(B)')
axs[2].axis('off')

plt.show()

③AᵀA を計算

# NumPyを使った行列積
M_r = np.dot(A_rT, A_r)  # または M_r = A_rT @ A_r
M_g = np.dot(A_gT, A_g)  # または M_g = A_gT @ A_g
M_b = np.dot(A_bT, A_b)  # または M_b = A_bT @ A_b
print("行列 M の形状:", M_r.shape)

④固有値の計算 & 固有ベクトルの計算

# 固有値を計算（NumPy の関数を使用せずに実装するのは長くなるため、簡略化します）
def approximate_eigenvalues(matrix):
    # 対角成分を固有値の近似として使用
    eigenvals, eigenvecs = np.linalg.eig(matrix)
    return eigenvals, eigenvecs

# 固有値を計算
eigenvals_r, eigenvecs_r = approximate_eigenvalues(M_r)
eigenvals_g, eigenvecs_g = approximate_eigenvalues(M_g)
eigenvals_b, eigenvecs_b = approximate_eigenvalues(M_b)

# 固有値をプロット
fig, axs = plt.subplots(1, 3, figsize=(20, 5))

axs[0].plot(eigenvals_r)
axs[0].set_title('Eigenvalues(R)')

axs[1].plot(eigenvals_g)
axs[1].set_title('Eigenvalues(G)')

axs[2].plot(eigenvals_b)
axs[2].set_title('Eigenvalues(B)')

plt.show()

⑤固有ベクトルの表示

# 固有ベクトルを表示（上位2つを2Dベクトルとしてプロット）
num_vectors = 2  # プロットする固有ベクトルの数

# 2D平面上に固有ベクトルをプロット
plt.figure(figsize=(6, 6))
origin = [0, 0]  # ベクトルの始点

for i in range(num_vectors):
    eigenvec_r = eigenvecs_r[:, i]
    # 固有ベクトルを正規化
    eigenvec_norm = eigenvec_r / np.linalg.norm(eigenvec_r)
    # 上位2つの成分を使用（高次元なので2Dに射影）
    vec = eigenvec_norm[:2]
    plt.quiver(*origin, vec[0], vec[1], angles='xy', scale_units='xy', scale=1, label=f'Eigenvector {i+1}')

plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.title('Eigenvectors')
plt.legend()
plt.grid()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

⑥特異値の計算(√λ (λ=固有値))

# 特異値を計算（固有値の平方根）
singular_values_r = np.sqrt(np.abs(eigenvals_r))
singular_values_g = np.sqrt(np.abs(eigenvals_g))
singular_values_b = np.sqrt(np.abs(eigenvals_b))

# 特異値をプロット
fig, axs = plt.subplots(1, 3, figsize=(20, 5))

axs[0].plot(singular_values_r)
axs[0].set_title('Singularvalues(R)')

axs[1].plot(singular_values_g)
axs[1].set_title('Singularvalues(G)')

axs[2].plot(singular_values_b)
axs[2].set_title('Singularvalues(B)')

plt.show()

⑦行列 V の構築

# 行列 V を構築（固有ベクトルを列として持つ）
V_r = eigenvecs_r
V_g = eigenvecs_g
V_b = eigenvecs_b

# 行列 V を画像として表示
fig, axs = plt.subplots(1, 3, figsize=(20, 5))

#axs[0].imshow(V_r)
axs[0].imshow(V_r.real)
axs[0].set_title('Matrix V (Eigenvectors)(R)')

#axs[1].imshow(V_g)
axs[1].imshow(V_g.real)
axs[1].set_title('Matrix V (Eigenvectors)(G)')

#axs[2].imshow(V_b)
axs[2].imshow(V_b.real)
axs[2].set_title('Matrix V (Eigenvectors)(B)')

plt.show()

⑧行列 U の計算(A=UΣVᵀ -> U=AVΣ⁻¹またはAAᵀの固有ベクトルで構築)

# Σ^(-1) を計算（特異値の逆数からなる対角行列）
Sigma_inv_r = np.diag([1/s if s > 1e-10 else 0 for s in singular_values_r])
Sigma_inv_g = np.diag([1/s if s > 1e-10 else 0 for s in singular_values_g])
Sigma_inv_b = np.diag([1/s if s > 1e-10 else 0 for s in singular_values_b])

# U = A * V * Σ^(-1) を計算
U_r = np.dot(A_r, np.dot(V_r, Sigma_inv_r))
U_g = np.dot(A_g, np.dot(V_g, Sigma_inv_g))
U_b = np.dot(A_b, np.dot(V_b, Sigma_inv_b))

# U の表示（正方形に切り取る）
min_size = min(U_r.shape[0], U_r.shape[1])
U_r_square = U_r[:min_size, :min_size]
U_g_square = U_g[:min_size, :min_size]
U_b_square = U_b[:min_size, :min_size]

# 行列 U を画像として表示
fig, axs = plt.subplots(1, 3, figsize=(20, 5))

axs[0].imshow(U_r_square.real)
axs[0].set_title('Matrix U (Left Singular Vectors)(R)')

axs[1].imshow(U_g_square.real)
axs[1].set_title('Matrix U (Left Singular Vectors)(G)')

axs[2].imshow(U_b_square.real)
axs[2].set_title('Matrix U (Left Singular Vectors)(B)')

plt.show()

⑨行列 A の再構築(A=UΣVᵀ)

# Σ（特異値からなる対角行列）
Sigma_r = np.diag(singular_values_r)
Sigma_g = np.diag(singular_values_g)
Sigma_b = np.diag(singular_values_b)

# V の転置を計算
V_rT = V_r.T
V_gT = V_g.T
V_bT = V_b.T

# U * Σ * V^T を用いて A を再構築
USigma_r = np.dot(U_r, Sigma_r)  
A_reconstructed_r = np.dot(USigma_r, V_rT)
USigma_g = np.dot(U_g, Sigma_g)  
A_reconstructed_g = np.dot(USigma_g, V_gT)
USigma_b = np.dot(U_b, Sigma_b)  
A_reconstructed_b = np.dot(USigma_b, V_bT)

# 再構築された画像を結合
reconstructed_img = np.stack((A_reconstructed_r.astype(np.uint8),
                                A_reconstructed_g.astype(np.uint8),
                                A_reconstructed_b.astype(np.uint8)), axis=2)


# オリジナル画像と再構築画像の差分を計算
difference = np.abs(image - reconstructed_img)

# 再構築された画像を表示
fig, axs = plt.subplots(1, 3, figsize=(10, 5))

axs[0].imshow(image, cmap='gray')
axs[0].set_title('Original Image')
axs[0].axis('off')

axs[1].imshow(reconstructed_img, cmap='gray')
axs[1].set_title('Reconstructed Image')
axs[1].axis('off')

# 差分画像を表示
axs[2].imshow(difference, cmap='gray')
axs[2].set_title('Difference Image')
axs[2].axis('off')
plt.show()

⑩画像圧縮

# ランクを削減する（k = 保持する特異値の数）
# ランクは特異値であるSの値であり、ランク1は1行目、ランク2は2行目の特異値を表しています
# SVDを用いて画像を圧縮する方法はランクを省いて行列を小さくする方法を使っているので、
# 一番大きい値を持つ1行目のランクからn行目までのランクを使う方式で圧縮が行われる
k_values = [10]  # 確認したい特異値の数(複数設定して各々ランク数に対する結果も得られる)

# プロットの準備
fig, axs = plt.subplots(len(k_values), 2, figsize=(12, len(k_values) * 6))

for idx, k in enumerate(k_values):
    # 特異値と対応する行列を削減
    Sigma_r_k = np.diag(singular_values_r[:k])
    Sigma_g_k = np.diag(singular_values_g[:k])
    Sigma_b_k = np.diag(singular_values_b[:k])

    U_r_k = U_r[:, :k]
    U_g_k = U_g[:, :k]
    U_b_k = U_b[:, :k]

    V_rT_k = V_rT[:k, :]
    V_gT_k = V_gT[:k, :]
    V_bT_k = V_bT[:k, :]

    # 削減されたランクで画像を再構築
    A_reconstructed_r_k = np.dot(np.dot(U_r_k, Sigma_r_k), V_rT_k)
    A_reconstructed_g_k = np.dot(np.dot(U_g_k, Sigma_g_k), V_gT_k)
    A_reconstructed_b_k = np.dot(np.dot(U_b_k, Sigma_b_k), V_bT_k)

    # 再構築された画像を結合
    reconstructed_img_k = np.stack((A_reconstructed_r_k,
                                     A_reconstructed_g_k,
                                     A_reconstructed_b_k), axis=2)

    # 値をクリップしてキャスト
    reconstructed_img_k = np.clip(reconstructed_img_k, 0, 255).astype(np.uint8)

    # 差分画像の計算
    diff_img = np.abs(image.astype(np.float32) - reconstructed_img_k.astype(np.float32)).astype(np.uint8)

    # プロット
    axs[idx, 0].imshow(reconstructed_img_k)
    axs[idx, 0].set_title(f'Reconstructed Image (Rank-{k})')
    axs[idx, 0].axis('off')

    axs[idx, 1].imshow(diff_img)
    axs[idx, 1].set_title(f'Difference (Rank-{k})')
    axs[idx, 1].axis('off')

plt.tight_layout()
plt.show()

結果

SVD圧縮, Jpeg, WebPを使って画像の容量を全体の約5%(350Kb)ぐらいに減らした結果、SVD圧縮の結果得られた３つの行列(U,S,V)の総要素数は92,151個で、SVDを用いて圧縮した行列の容量は4Bytes x [ 画像幅(2184) * ランク数(27) + ランク数(27) + 画像高さ(1228) * ランク数(27) ] = 368,604Bytes(約368KB)でした。VD圧縮で使用された特異値の数は9個でした。（元の画像の容量は 3Bytes x 画像幅(2184) × 画像高さ(1228) = 8,045,856 Bytes）

そして、SVD圧縮で減らした容量に合わせて、JpegとWebPも損失圧縮し、Jpegの容量は364KB、WebPの容量は366KBにして画像の品質を比較した結果SVDの品質はノイズが多く画像の重要な情報が多く損失されていることを確認できました。

そして元の画像に対するSVD圧縮、Jpeg、WebPのPSNR、SSIM、MSEを計算した結果、PSNRはWebP>Jpeg>SVD、SSIMはJpeg>WebP>SVD、MSEはJpeg<WebP<SVDでした。SSIMはJpegとWebPが大きな差がなかったため、圧縮の品質としてはWebP>Jpeg>SVDの順でした。

Jpegの容量は364KB、WebPの容量は366KBの画質に近似できる品質を得るためにはSVDを適用した情報量を削減する際に380個のランクを使う必要がありました、この場合SVDを用いて圧縮した画像のPSNR、SSIM、MSEはJpegとWebPの評価値と近似できますが、容量は4,340KBになり、比較的に多い容量が必要でした。

結論

SVDを用いた画像圧縮の有用性としては、画像の特徴を基にして圧縮ができる点と、入力画像の特異値の数を調整することで圧縮率と画質を柔軟に制御できる点は、大きな利点です。

その反面、実用性の観点からはいくつかの課題が存在します。SVDを用いた画像圧縮はJpegやWebPなどの一般的な画像圧縮アルゴリズムと比較すると、圧縮後の画質が劣化しやすく、画像圧縮という特定の用途において競争力が低いといえます。

総じて、SVDによる画像圧縮は特徴を基にした圧縮やその圧縮率を柔軟に調整できるのは利点であるが、一般的な画像圧縮用途での実用性は限定的であると評価されます。

以上、ありがとうございました！