使用核心 API 進行矩陣近似

在 TensorFlow.org 上檢視 在 Google Colab 中執行 在 GitHub 上檢視原始碼 下載筆記本

簡介

這個筆記本使用 TensorFlow Core 低階 API 來展示 TensorFlow 作為高效能科學運算平台的強大功能。請造訪「Core API 總覽」以深入瞭解 TensorFlow Core 及其預期用途。

本教學課程探討奇異值分解 (SVD) 技術及其在低秩近似問題中的應用。SVD 用於分解實數或複數矩陣，並在資料科學中有各種用途，例如圖片壓縮。本教學課程的圖片來自 Google Brain 的 Imagen 專案。

svd_intro

設定

import matplotlib
from matplotlib.image import imread
from matplotlib import pyplot as plt
import requests
# Preset Matplotlib figure sizes.
matplotlib.rcParams['figure.figsize'] = [16, 9]

import tensorflow as tf
print(tf.__version__)

SVD 基本原理

矩陣 \({\mathrm{A} }\) 的奇異值分解由以下因式分解決定

\[{\mathrm{A} } = {\mathrm{U} } \Sigma {\mathrm{V} }^T\]

其中

  • \(\underset{m \times n}{\mathrm{A} }\)：輸入矩陣，其中 \(m \geq n\)
  • \(\underset{m \times n}{\mathrm{U} }\)：正交矩陣，\({\mathrm{U} }^T{\mathrm{U} } = {\mathrm{I} }\)，其中每個資料欄 \(u_i\) 都表示 \({\mathrm{A} }\) 的左奇異向量
  • \(\underset{n \times n}{\Sigma}\)：對角矩陣，其中每個對角線項目 \(\sigma_i\) 都表示 \({\mathrm{A} }\) 的奇異值
  • \(\underset{n \times n}{ {\mathrm{V} }^T}\)：正交矩陣，\({\mathrm{V} }^T{\mathrm{V} } = {\mathrm{I} }\)，其中每個資料列 \(v_i\) 都表示 \({\mathrm{A} }\) 的右奇異向量

當 \(m < n\) 時，\({\mathrm{U} }\) 和 \(\Sigma\) 的維度皆為 \((m \times m)\)，且 \({\mathrm{V} }^T\) 的維度為 \((m \times n)\)。

svd_full

TensorFlow 的線性代數套件具有函式 tf.linalg.svd，可用於計算一個或多個矩陣的奇異值分解。首先定義一個簡單的矩陣並計算其 SVD 因式分解。

A = tf.random.uniform(shape=[40,30])
# Compute the SVD factorization
s, U, V = tf.linalg.svd(A)
# Define Sigma and V Transpose
S = tf.linalg.diag(s)
V_T = tf.transpose(V)
# Reconstruct the original matrix
A_svd = U@S@V_T
# Visualize 
plt.bar(range(len(s)), s);
plt.xlabel("Singular value rank")
plt.ylabel("Singular value")
plt.title("Bar graph of singular values");

tf.einsum 函式可用於直接從 tf.linalg.svd 的輸出計算矩陣重建。

A_svd = tf.einsum('s,us,vs -> uv',s,U,V)
print('\nReconstructed Matrix, A_svd', A_svd)

使用 SVD 進行低秩近似

矩陣 \({\mathrm{A} }\) 的秩由其資料欄所跨向量空間的維度決定。SVD 可用於以較低的秩近似矩陣，最終減少儲存矩陣所表示資訊所需的資料維度。

SVD 中的 \({\mathrm{A} }\) 秩 r 近似值由以下公式定義

\[{\mathrm{A_r} } = {\mathrm{U_r} } \Sigma_r {\mathrm{V_r} }^T\]

其中

  • \(\underset{m \times r}{\mathrm{U_r} }\)：由 \({\mathrm{U} }\) 的前 \(r\) 個資料欄組成的矩陣
  • \(\underset{r \times r}{\Sigma_r}\)：由 \(\Sigma\) 中的前 \(r\) 個奇異值組成的對角矩陣
  • \(\underset{r \times n}{\mathrm{V_r} }^T\): 由 \({\mathrm{V} }^T\) 的前 \(r\) 個資料列組成的矩陣

svd_approx

首先編寫一個函式，以計算給定矩陣的秩 r 近似值。此低秩近似程序用於圖片壓縮；因此，計算每個近似值的實際資料大小也很有幫助。為簡化起見，假設秩 r 近似矩陣的資料大小等於計算近似值所需的元素總數。接著，編寫一個函式以視覺化原始矩陣 \(\mathrm{A}\)、其秩 r 近似值 \(\mathrm{A}_r\) 和誤差矩陣 \(|\mathrm{A} - \mathrm{A}_r|\)。

def rank_r_approx(s, U, V, r, verbose=False):
  # Compute the matrices necessary for a rank-r approximation
  s_r, U_r, V_r = s[..., :r], U[..., :, :r], V[..., :, :r] # ... implies any number of extra batch axes
  # Compute the low-rank approximation and its size
  A_r = tf.einsum('...s,...us,...vs->...uv',s_r,U_r,V_r)
  A_r_size = tf.size(U_r) + tf.size(s_r) + tf.size(V_r)
  if verbose:
    print(f"Approximation Size: {A_r_size}")
  return A_r, A_r_size

def viz_approx(A, A_r):
  # Plot A, A_r, and A - A_r
  vmin, vmax = 0, tf.reduce_max(A)
  fig, ax = plt.subplots(1,3)
  mats = [A, A_r, abs(A - A_r)]
  titles = ['Original A', 'Approximated A_r', 'Error |A - A_r|']
  for i, (mat, title) in enumerate(zip(mats, titles)):
    ax[i].pcolormesh(mat, vmin=vmin, vmax=vmax)
    ax[i].set_title(title)
    ax[i].axis('off')

print(f"Original Size of A: {tf.size(A)}")
s, U, V = tf.linalg.svd(A)

# Rank-15 approximation
A_15, A_15_size = rank_r_approx(s, U, V, 15, verbose = True)
viz_approx(A, A_15)

# Rank-3 approximation
A_3, A_3_size = rank_r_approx(s, U, V, 3, verbose = True)
viz_approx(A, A_3)

如同預期，使用較低的秩會導致近似值精確度降低。不過，在實際情況中，這些低秩近似值的品質通常已足夠。另請注意，使用 SVD 進行低秩近似的主要目標是減少資料的維度，而非減少資料本身的磁碟空間。然而，隨著輸入矩陣的維度變得更高，許多低秩近似值最終也會受益於縮減的資料大小。這種縮減效益是此程序適用於圖片壓縮問題的原因。

圖片載入

以下圖片可在 Imagen 首頁上找到。Imagen 是 Google Research 的 Brain 團隊開發的文字轉圖片擴散模型。人工智慧根據提示「在時代廣場騎自行車的柯基犬照片。牠戴著太陽眼鏡和沙灘帽」建立了這張圖片。是不是很酷！您也可以將下方的網址變更為任何 .jpg 連結，以載入自選的自訂圖片。

首先讀取並視覺化圖片。讀取 JPEG 檔案後，Matplotlib 會輸出形狀為 \((m \times n \times 3)\) 的矩陣 \({\mathrm{I} }\)，代表具有紅色、綠色和藍色這 3 個顏色通道的 2 維圖片。

img_link = "https://imagen.research.google/main_gallery_images/a-photo-of-a-corgi-dog-riding-a-bike-in-times-square.jpg"
img_path = requests.get(img_link, stream=True).raw
I = imread(img_path, 0)
print("Input Image Shape:", I.shape)

def show_img(I):
  # Display the image in matplotlib
  img = plt.imshow(I)
  plt.axis('off')
  return

show_img(I)

圖片壓縮演算法

現在，使用 SVD 計算範例圖片的低秩近似值。請注意，圖片的形狀為 \((1024 \times 1024 \times 3)\)，且理論 SVD 僅適用於 2 維矩陣。這表示範例圖片必須分批處理成 3 個大小相等的矩陣，分別對應至 3 個顏色通道。您可以將矩陣轉置為形狀 \((3 \times 1024 \times 1024)\) 來達成此目的。為了清楚視覺化近似誤差，請將圖片的 RGB 值從 \([0,255]\) 重新調整為 \([0,1]\)。請記得先將近似值裁剪至此區間內，再進行視覺化。tf.clip_by_value 函式對此很有幫助。

def compress_image(I, r, verbose=False):
  # Compress an image with the SVD given a rank 
  I_size = tf.size(I)
  print(f"Original size of image: {I_size}")
  # Compute SVD of image
  I = tf.convert_to_tensor(I)/255
  I_batched = tf.transpose(I, [2, 0, 1]) # einops.rearrange(I, 'h w c -> c h w')
  s, U, V = tf.linalg.svd(I_batched)
  # Compute low-rank approximation of image across each RGB channel
  I_r, I_r_size = rank_r_approx(s, U, V, r)
  I_r = tf.transpose(I_r, [1, 2, 0]) # einops.rearrange(I_r, 'c h w -> h w c')
  I_r_prop = (I_r_size / I_size)
  if verbose:
    # Display compressed image and attributes
    print(f"Number of singular values used in compression: {r}")
    print(f"Compressed image size: {I_r_size}")
    print(f"Proportion of original size: {I_r_prop:.3f}")
    ax_1 = plt.subplot(1,2,1)
    show_img(tf.clip_by_value(I_r,0.,1.))
    ax_1.set_title("Approximated image")
    ax_2 = plt.subplot(1,2,2)
    show_img(tf.clip_by_value(0.5+abs(I-I_r),0.,1.))
    ax_2.set_title("Error")
  return I_r, I_r_prop

現在，計算以下秩的秩 r 近似值：100、50、10

I_100, I_100_prop = compress_image(I, 100, verbose=True)

I_50, I_50_prop = compress_image(I, 50, verbose=True)

I_10, I_10_prop = compress_image(I, 10, verbose=True)

評估近似值

有多種有趣的方法可以衡量有效性，並更妥善掌控矩陣近似值。

壓縮率與秩的比較

針對上述每個近似值，觀察資料大小如何隨秩變化。

plt.figure(figsize=(11,6))
plt.plot([100, 50, 10], [I_100_prop, I_50_prop, I_10_prop])
plt.xlabel("Rank")
plt.ylabel("Proportion of original image size")
plt.title("Compression factor vs rank");

根據此圖表，近似圖片的壓縮率與其秩之間存在線性關係。為了進一步探索這一點，請注意近似矩陣 \({\mathrm{A} }_r\) 的資料大小定義為計算所需的元素總數。以下方程式可用於找出壓縮率與秩之間的關係

\[x = (m \times r) + r + (r \times n) = r \times (m + n + 1)\]

\[c = \large \frac{x}{y} = \frac{r \times (m + n + 1)}{m \times n}\]

其中

  • \(x\)：\({\mathrm{A_r} }\) 的大小
  • \(y\)：\({\mathrm{A} }\) 的大小
  • \(c = \frac{x}{y}\)：壓縮率
  • \(r\)：近似值的秩
  • \(m\) 和 \(n\)：\({\mathrm{A} }\) 的資料列和資料欄維度

為了找出將圖片壓縮到所需比例 \(c\) 所需的秩 \(r\)，可以重新排列上述方程式以解出 \(r\)

\[r = ⌊{\large\frac{c \times m \times n}{m + n + 1} }⌋\]

請注意，此公式與顏色通道維度無關，因為每個 RGB 近似值彼此之間互不影響。現在，編寫一個函式，以在給定所需壓縮率的情況下壓縮輸入圖片。

def compress_image_with_factor(I, compression_factor, verbose=False):
  # Returns a compressed image based on a desired compression factor
  m,n,o = I.shape
  r = int((compression_factor * m * n)/(m + n + 1))
  I_r, I_r_prop = compress_image(I, r, verbose=verbose)
  return I_r

將圖片壓縮至原始大小的 15%。

compression_factor = 0.15
I_r_img = compress_image_with_factor(I, compression_factor, verbose=True)

奇異值的累積總和

奇異值的累積總和可有效指出秩 r 近似值所擷取的能量。視覺化範例圖片中 RGB 平均的奇異值累積比例。tf.cumsum 函式對此很有幫助。

def viz_energy(I):
  # Visualize the energy captured based on rank
  # Computing SVD
  I = tf.convert_to_tensor(I)/255
  I_batched = tf.transpose(I, [2, 0, 1]) 
  s, U, V = tf.linalg.svd(I_batched)
  # Plotting average proportion across RGB channels 
  props_rgb = tf.map_fn(lambda x: tf.cumsum(x)/tf.reduce_sum(x), s)
  props_rgb_mean = tf.reduce_mean(props_rgb, axis=0)
  plt.figure(figsize=(11,6))
  plt.plot(range(len(I)), props_rgb_mean, color='k')
  plt.xlabel("Rank / singular value number")
  plt.ylabel("Cumulative proportion of singular values")
  plt.title("RGB-averaged proportion of energy captured by the first 'r' singular values")

viz_energy(I)

看起來這張圖片中超過 90% 的能量都擷取在前 100 個奇異值內。現在，編寫一個函式，以在給定所需能量保留率的情況下壓縮輸入圖片。

def compress_image_with_energy(I, energy_factor, verbose=False):
  # Returns a compressed image based on a desired energy factor
  # Computing SVD
  I_rescaled = tf.convert_to_tensor(I)/255
  I_batched = tf.transpose(I_rescaled, [2, 0, 1]) 
  s, U, V = tf.linalg.svd(I_batched)
  # Extracting singular values
  props_rgb = tf.map_fn(lambda x: tf.cumsum(x)/tf.reduce_sum(x), s)
  props_rgb_mean = tf.reduce_mean(props_rgb, axis=0)
  # Find closest r that corresponds to the energy factor
  r = tf.argmin(tf.abs(props_rgb_mean - energy_factor)) + 1
  actual_ef = props_rgb_mean[r]
  I_r, I_r_prop = compress_image(I, r, verbose=verbose)
  print(f"Proportion of energy captured by the first {r} singular values: {actual_ef:.3f}")
  return I_r

將圖片壓縮至保留 75% 的能量。

energy_factor = 0.75
I_r_img = compress_image_with_energy(I, energy_factor, verbose=True)

誤差與奇異值

近似誤差與奇異值之間也存在有趣的關係。結果顯示，近似值的平方 Frobenius 範數等於遺漏的奇異值平方和

\[{||A - A_r||}^2 = \sum_{i=r+1}^{R}σ_i^2\]

使用本教學課程一開始的範例矩陣的秩 10 近似值，測試此關係。

s, U, V = tf.linalg.svd(A)
A_10, A_10_size = rank_r_approx(s, U, V, 10)
squared_norm = tf.norm(A - A_10)**2
s_squared_sum = tf.reduce_sum(s[10:]**2)
print(f"Squared Frobenius norm: {squared_norm:.3f}")
print(f"Sum of squared singular values left out: {s_squared_sum:.3f}")

結論

這個筆記本介紹了使用 TensorFlow 實作奇異值分解，並將其應用於編寫圖片壓縮演算法的流程。以下是一些可能有幫助的訣竅

  • TensorFlow Core API 可用於各種高效能科學運算用途。
  • 若要進一步瞭解 TensorFlow 的線性代數功能，請造訪 linalg 模組的文件。
  • SVD 也可用於建構推薦系統。

如需更多使用 TensorFlow Core API 的範例，請查看指南。如果您想進一步瞭解如何載入和準備資料，請參閱圖片資料載入或CSV 資料載入教學課程。