pylops

PyLops - 線形演算子ライブラリ

クイックリファレンス

import numpy as np
import pylops

# Create operator and apply forward/adjoint
A = pylops.FirstDerivative(n=100, dtype='float64')
y = A @ x        # Forward: y = A @ x
x_adj = A.H @ y  # Adjoint: x = A.H @ y
x_est = A / y    # Solve inverse problem

主なクラス

クラス	用途
LinearOperator	すべての演算子の基底クラス
VStack/HStack	演算子を垂直／水平方向にスタック
BlockDiag	ブロック対角演算子の合成

基本的な操作

基本演算子

# Diagonal operator
D = pylops.Diagonal(np.array([1., 2., 3.]))
y = D @ x; x_adj = D.H @ y

# Derivatives
D1 = pylops.FirstDerivative(n, dtype='float64')
D2 = pylops.SecondDerivative(n, dtype='float64')
G = pylops.Gradient(dims=(64, 64), dtype='float64')

畳み込み

wavelet = np.sin(np.linspace(0, 2*np.pi, 21)) * np.hanning(21)
C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=10)
y = C @ x      # Convolve
x_adj = C.H @ y  # Correlation (adjoint)

演算子の合成とスタック

# Chain: y = C @ B @ A @ x
composed = pylops.Smoothing1D(5, n) @ pylops.FirstDerivative(n) @ pylops.Identity(n)

# Stack operators
V = pylops.VStack([A, B])     # Vertical: (2n, n)
H = pylops.HStack([A, B])     # Horizontal: (n, 2n)
BD = pylops.BlockDiag([A, B]) # Block diagonal: (2n, 2n)

逆問題を解く

# Simple least squares
x_est = A / y

# Normal equations
x_est = pylops.optimization.leastsquares.NormalEquationsInversion(A, None, y)

# Regularized inversion with smoothness
Reg = pylops.SecondDerivative(n)
x_est = pylops.optimization.leastsquares.RegularizedInversion(
    A, [Reg], y, epsRs=[0.1]
)

反復ソルバー

x_lsqr = pylops.optimization.solver.lsqr(A, y, iter_lim=100)[0]
x_cgls = pylops.optimization.solver.cgls(A, y, niter=100)[0]

スパース性を促進する逆変換

x_l1 = pylops.optimization.sparsity.fista(A, y, niter=100, eps=0.1)[0]

随伴演算子の検証（ドットテスト）

A = pylops.FirstDerivative(100)
pylops.utils.dottest(A, 100, 100, verb=True)  # Dot test passed!

よくあるパターン

地震波デコンボリューション

C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=len(wavelet)//2)
seismic = C @ reflectivity

# Deconvolve (regularized)
Reg = pylops.SecondDerivative(n)
reflectivity_est = pylops.optimization.leastsquares.RegularizedInversion(
    C, [Reg], seismic, epsRs=[0.01]
)

TV による画像のノイズ除去

ny, nx = image.shape
G = pylops.Gradient(dims=(ny, nx))
x_tv = pylops.optimization.sparsity.splitbregman(
    pylops.Identity(ny*nx), G, image.ravel(),
    niter_inner=5, niter_outer=10, mu=1.0, epsRL1s=[0.1]
)[0].reshape(ny, nx)

代替手段との使い分け

シナリオ	推奨
大規模な逆問題向けの行列を明示しない線形演算子	PyLops - 専用設計でメモリ効率に優れています
構造が既知のスパース行列演算	scipy.sparse - 標準的で、ドキュメントが充実しています
単純な畳み込み／デコンボリューション	PyLops - Convolve1D による使いやすい API
小規模な問題向けのカスタム演算子	Custom NumPy/SciPy - 追加の依存関係は不要です
GPU アクセラレーション対応の線形代数	PyLops - CuPy 配列を渡すだけで GPU を自動的に利用できます
地震波デコンボリューションやイメージング演算子	PyLops - 信号処理演算子の豊富なライブラリ

PyLops を選ぶ場合：明示的な行列を作成せずに、大規模な問題にも対応できる線形演算子が必要な場合です。演算子代数（@、VStack、BlockDiag）と組み込みソルバー（LSQR、FISTA、Split Bregman）により、逆問題のワークフローを簡潔に記述できます。

PyLops を避ける場合：明示的な行列で扱えるほど問題が小さい場合（NumPy/SciPy を使用してください）、または非線形演算子が必要な場合です（PyLops が扱うのは線形演算子のみです）。

よくあるワークフロー

正則化付き地震波デコンボリューション

• ウェーブレット配列を定義し、Convolve1D 演算子を作成する
• 地震波トレースデータを生成または読み込む
• ドットテストを実行し、演算子の随伴を検証する：pylops.utils.dottest()
• 正則化演算子（例：滑らかさを表す SecondDerivative）を設定する
• RegularizedInversion(C, [Reg], data, epsRs=[eps]) を実行する
• 推定した反射率を、真の値（利用可能な場合）と比較する
• epsRs パラメーターを調整する：値が大きいほど滑らかに、小さいほどシャープになります
• スパースな解には、代わりに pylops.optimization.sparsity.fista() を使用する

ヒント

1. 完全な行列を作成しない - メモリ効率を高めるため、.matvec() と .rmatvec() を使用します
2. 形状を確認する - 演算子には行列と同様に .shape 属性があります
3. 随伴を検証する - カスタム演算子では必ず pylops.utils.dottest() を実行します
4. 簡単なものから始める - 規模を拡大する前に、小規模な問題でテストします
5. GPU を使用する - CuPy 配列を渡すと、GPU アクセラレーションが自動的に有効になります

参考資料

• 演算子リファレンス - 利用可能な演算子の完全な一覧
• ソルバーリファレンス - ソルバーの手法と設定オプション

スクリプト

• scripts/deconvolution.py - 地震波デコンボリューションの例