Juliaで学ぶ線形代数：reflect!()を使った鏡面反射のプログラミング

2025-03-21

JuliaのLinearAlgebra.reflect!()関数は、与えられたベクトルを軸とした鏡面反射（または鏡像変換）を、対象のベクトルまたは行列に対して破壊的に実行する関数です。 "破壊的"というのは、元のオブジェクト自体が変更されるという意味です。

もう少し詳しく説明します。

鏡面反射とは

鏡面反射とは、ある超平面（2次元なら直線、3次元なら平面）を鏡として、その鏡に対して対称な像を作る操作のことです。 LinearAlgebra.reflect!()では、指定されたベクトルがこの超平面の法線ベクトル（超平面に垂直なベクトル）となります。

LinearAlgebra.reflect!() の使い方

reflect!(v, x)

x: 反射の対象となるベクトルまたは行列。
v: 反射の軸となる超平面の法線ベクトル。

具体例

using LinearAlgebra

# 反射の軸となる法線ベクトル
v = [1.0, 0.0]

# 反射されるベクトル
x = [2.0, 1.0]

# 反射を実行 (xが変更される)
reflect!(v, x)

println(x)  # 出力: [-2.0, 1.0]

# 行列の場合
A = [2.0 1.0; 3.0 4.0]
reflect!(v, A)
println(A) # 出力: [-2.0 1.0; -3.0 4.0]

上記の例では、v = [1.0, 0.0] がx軸を鏡とした反射を表しています。 x = [2.0, 1.0] は、x軸を鏡として反射されると [-2.0, 1.0] になります。行列Aの場合も同様に、各列ベクトルが反射されます。

v の大きさは反射の結果に影響しません。v の方向のみが重要です。
v はゼロベクトルであってはいけません。
reflect!() は破壊的な関数なので、元のxが変更されます。元のxを保持しておきたい場合は、copy(x) などでコピーを作成してから reflect!() を使用してください。

v (反射軸ベクトル) がゼロベクトル

対処法
v にゼロベクトル以外の値を設定してください。例えば、x軸方向の反射なら v = [1.0, 0.0]、y軸方向なら v = [0.0, 1.0] のように設定します。
理由
ゼロベクトルは方向を持たないため、鏡面反射を定義できません。
間違い
v = [0.0, 0.0] のように、反射軸ベクトルがゼロベクトルになっている。

x (反射対象) の型が不適切

対処法
x をベクトルまたは行列に変換してください。必要であれば、convert() 関数などを使用します。
理由
reflect!() は、ベクトルまたは行列に対してのみ定義されています。
間違い
x がベクトルや行列以外の型（例えば、整数や文字列）である。

x の次元と v の次元が一致しない (ベクトルに対して)

対処法
x と v の次元を一致させてください。例えば、x が2次元ベクトルなら、v も2次元ベクトルである必要があります。
理由
反射は、x と同じ次元空間内で行われる必要があります。
間違い
x が2次元ベクトルなのに、v が3次元ベクトルである。

行列に対する反射の解釈

対処法
行列の反射がどのように行われるかを理解しておくことが重要です。
理由
reflect!(v, A) は、行列 A の各列ベクトルを個別に v を軸として反射します。
間違い
行列 A に対して reflect!(v, A) を実行したとき、行列全体がひとつのオブジェクトとして反射されると誤解している。

破壊的変更の影響

対処法
x のコピーを作成してから reflect!() を使用してください。例えば、x_copy = copy(x) として、reflect!(v, x_copy) を実行します。
理由
reflect!() は破壊的な関数なので、x 自体が変更されます。
間違い
元の x の値を保持しておきたいのに、reflect!() を使ってしまい、元の値が失われた。

数値誤差

対処法
許容範囲内の誤差であれば、そのまま使用するか、必要に応じて近似的な計算を行います。
理由
浮動小数点数演算には、わずかな誤差が常に存在します。
間違い
浮動小数点数の計算誤差により、期待通りの結果が得られない。

Juliaのバージョン

対処法
最新の安定版Juliaを使用することを推奨します。
理由
Juliaのバージョンによって、関数やライブラリの仕様が変更されることがあります。
間違い
古いバージョンのJuliaを使用しているため、reflect!() が存在しない、または挙動が異なる。

エラーメッセージをよく読む
Juliaが表示するエラーメッセージは、問題解決のヒントになります。
typeof() で型を確認
変数の型が意図したものと一致しているか確認します。
println() で値を確認
変数の値が期待通りになっているか確認します。
ドキュメントを参照
Juliaのドキュメントで、reflect!() の使い方や注意点を確認します。
最小限のコードで試す
問題を切り分けるために、できるだけ短いコードで問題を再現させます。

2次元ベクトルの反射 (x軸を軸とする)

using LinearAlgebra

# 反射軸 (x軸)
v = [1.0, 0.0]

# 反射されるベクトル
x = [2.0, 1.0]

# 反射を実行
reflect!(v, x)

println(x)  # 出力: [-2.0, 1.0]

この例では、v = [1.0, 0.0] がx軸を表す法線ベクトルです。x = [2.0, 1.0] は、x軸を鏡として反射され、[-2.0, 1.0] になります。

2次元ベクトルの反射 (任意の軸)

using LinearAlgebra

# 反射軸 (任意)
v = [cos(pi/4), sin(pi/4)] # 45度傾いた軸

# 反射されるベクトル
x = [1.0, 0.0]

# 反射を実行
reflect!(v, x)

println(x)

この例では、v が45度傾いた軸を表す法線ベクトルです。 cos(pi/4)とsin(pi/4)で、単位ベクトルを生成しています。

3次元ベクトルの反射

using LinearAlgebra

# 反射軸 (z軸)
v = [0.0, 0.0, 1.0]

# 反射されるベクトル
x = [1.0, 2.0, 3.0]

# 反射を実行
reflect!(v, x)

println(x)  # 出力: [1.0, 2.0, -3.0]

この例では、z軸を鏡面として反射が行われます。

行列の反射

using LinearAlgebra

# 反射軸 (x軸)
v = [1.0, 0.0]

# 反射される行列
A = [1.0 2.0; 3.0 4.0]

# 反射を実行
reflect!(v, A)

println(A) # 出力: [-1.0 -2.0; 3.0 4.0]

この例では、行列 A の各列ベクトルがx軸を鏡として反射されます。

破壊的変更を避ける

using LinearAlgebra

# 反射軸 (x軸)
v = [1.0, 0.0]

# 反射されるベクトル (コピーを作成)
x = [2.0, 1.0]
x_copy = copy(x)

# コピーに対して反射を実行
reflect!(v, x_copy)

println(x)      # 出力: [2.0, 1.0] (元のxは変更されない)
println(x_copy) # 出力: [-2.0, 1.0] (コピーは変更される)

この例では、copy() 関数を使って x のコピーを作成し、コピーに対して reflect!() を実行しています。これにより、元の x の値が変更されるのを防ぐことができます。

using LinearAlgebra

function reflect_over_hyperplane!(x, n)
    # n: 超平面の法線ベクトル (単位ベクトルである必要はない)
    n = n / norm(n) # 法線ベクトルを単位ベクトルに正規化
    reflect!(n, x)
end

x = [1.0, 2.0, 3.0]
n = [1.0, 1.0, 0.0] # 超平面の法線ベクトル
reflect_over_hyperplane!(x, n)
println(x)

手動での計算

最も基本的な方法は、反射の計算を直接行うことです。これは、reflect!() の内部で行われている計算を自分で実装するものです。

function reflect_manual!(x, v)
  v = v / norm(v) # vを単位ベクトルに正規化 (重要!)
  dot_xv = dot(x, v)
  x[:] = x - 2 * dot_xv * v # xを破壊的に更新
end

x = [1.0, 2.0]
v = [1.0, 0.0]
reflect_manual!(x, v)
println(x) # 出力: [-1.0, 2.0]

#行列の場合
A = [1.0 2.0; 3.0 4.0]
for i in 1:size(A,2)
    reflect_manual!(view(A,:,i), v) #各列ベクトルに対して処理
end
println(A) # 出力: [-1.0 -2.0; 3.0 4.0]

この方法の利点は、反射の計算を完全に理解できることです。欠点は、コードが少し長くなることと、vを自分で正規化する必要があることです。reflect!()は内部で正規化を行っています。行列の場合、各列ベクトルに対して処理を行う必要があります。viewを使って元の行列のviewを渡すことで、メモリ効率の良い処理が可能です。

変換行列の利用

反射変換は、変換行列を使って表現することができます。この行列をベクトルや行列に掛けることで、反射の結果を得ることができます。

function reflection_matrix(v)
  v = v / norm(v)
  I = Matrix{Float64}(I, length(v), length(v)) # 単位行列
  return I - 2 * v * v' # 'は転置
end

function reflect_matrix!(x, v)
  R = reflection_matrix(v)
  x[:] = R * x  # xを破壊的に更新
end


x = [1.0, 2.0]
v = [1.0, 0.0]
reflect_matrix!(x, v)
println(x) # 出力: [-1.0, 2.0]

A = [1.0 2.0; 3.0 4.0]
for i in 1:size(A,2)
    reflect_matrix!(view(A,:,i), v)
end
println(A) # 出力: [-1.0 -2.0; 3.0 4.0]

この方法の利点は、変換行列を事前に計算しておけば、複数のベクトルや行列に対して効率的に反射を適用できることです。また、変換行列を他の変換と組み合わせることも容易になります。欠点は、reflect!()に比べてコードが少し長くなることです。

Householder変換の利用

Householder変換は、反射変換を一般化したものです。reflect!()もHouseholder変換の一種とみなすことができます。より複雑な反射操作を行う場合、Householder変換を直接使うことが有効です。JuliaのLinearAlgebraパッケージには、Householder変換に関連する関数が用意されています。

using LinearAlgebra

function reflect_householder!(x, v)
    H = householder_reflection(v)
    x[:] = H * x
end

x = [1.0, 2.0]
v = [1.0, 0.0]
reflect_householder!(x, v)
println(x)

A = [1.0 2.0; 3.0 4.0]
for i in 1:size(A,2)
    reflect_householder!(view(A,:,i), v)
end
println(A)

この例では、householder_reflection関数を使ってHouseholder反射オブジェクトを作成し、それをベクトルxに適用しています。

より複雑な反射操作を行う場合は、Householder変換を検討してください。
複数のベクトルや行列に対して同じ反射を適用する場合は、変換行列を使うと効率的です。
反射の計算を理解したい場合は、手動での計算が良いでしょう。
単純な反射であれば、reflect!() が最も手軽です。