Skip to content

Implement a robust tracking method based on Ogita-Aishima iteration #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
t-uda merged 2 commits into from
Sep 10, 2025
Merged

Implement a robust tracking method based on Ogita-Aishima iteration #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f382142
set Ogita-Aishima Iteration as default
Sep 9, 2025
3ebebda
add example large_n = 500
Sep 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
417 changes: 187 additions & 230 deletions notebooks/large_matrix_tracking_demo.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

74 changes: 53 additions & 21 deletions src/eigenpairflow/tracking.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,34 @@
from types import SimpleNamespace

from .ode import symmetric_ode_derivative
from .correction import correct_trajectory
from .correction import correct_trajectory, ogita_aishima_refinement
from .types import EigenTrackingResults


def _track_ogita_aishima(A_func, t_eval, *, max_iter=10, tol=1.0e-12, rho=1.0):
A_t = A_func(t_eval[0])
eigvals, eigvecs = scipy.linalg.eigh(A_t)
sort_indices = np.argsort(eigvals)
eigvals = eigvals[sort_indices]
eigvecs = eigvecs[:, sort_indices]
eigvals_list = [np.diag(eigvals)]
eigvecs_list = [eigvecs]
y0 = np.concatenate([eigvecs.flatten(), eigvals])
y = np.zeros((len(t_eval), len(y0)))
y[0] = y0
for i, t in enumerate(t_eval[1:], start=1):
A_t = A_func(t)
eigvecs, D = ogita_aishima_refinement(
A_t, eigvecs, max_iter=max_iter, tol=tol, rho=rho
)
eigvals_list.append(D)
eigvecs_list.append(eigvecs)
y[i] = np.concatenate([eigvecs.flatten(), np.diag(D)])

sol = SimpleNamespace(t=t_eval, y=y, success=True, message="Ogita-Aishima")
return eigvecs_list, eigvals_list, sol


def _track_symmetric_eigh(
A_func,
dA_func,
Expand Down Expand Up @@ -59,7 +83,9 @@ def _track_symmetric_eigh(
dt = t_values[i + 1] - t_i
dydt = symmetric_ode_derivative(t_i, y[:, i], n, dA_func)
y[:, i + 1] = y[:, i] + dt * dydt
sol = SimpleNamespace(t=t_values, y=y, success=True, message="Euler integration")
sol = SimpleNamespace(
t=t_values, y=y, success=True, message="Euler integration"
)
else:
sol = scipy.integrate.solve_ivp(
symmetric_ode_derivative,
Expand Down Expand Up @@ -88,8 +114,8 @@ def eigenpairtrack(
t_eval,
matrix_type="symmetric",
method="eigh",
correction_method="ogita_aishima",
solver_method="Euler",
solver_method="ogita_aishima",
correction_method=None,
Comment on lines 115 to +118

@chatgpt-codex-connector chatgpt-codex-connector Bot Sep 10, 2025

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[P1] Restore positional argument compatibility

The signature now lists solver_method before correction_method, whereas the previous release exposed the opposite order. Any existing call that passed these optional arguments positionally (e.g. eigenpairtrack(..., "matching", "RK45")) will now treat the correction value as the solver method and vice versa, typically resulting in invalid method names being forwarded to solve_ivp or correct_trajectory. Unless positional arguments are explicitly unsupported, this is a backwards‑incompatible change that will break working client code. Consider keeping the original parameter order or enforcing keyword‑only parameters to avoid silent API breakage.

Useful? React with 👍 / 👎.

rtol=1e-13,
atol=1e-12,
dense_output=False,
Expand All @@ -104,13 +130,16 @@ def eigenpairtrack(
dA_func (callable): 時刻 t に対して微分 dA/dt を返す関数。
t_span (tuple): 追跡する時間区間 (t_start, t_end)。
t_eval (np.ndarray): 結果を保存する時刻の配列。
matrix_type (str): 行列の種類。現在は ``"symmetric"`` のみ対応。
method (str): 固有分解の手法。現在は ``"eigh"`` のみ対応。
correction_method (str or None): 補正手法。 ``'ogita_aishima'`` がデフォルト。
solver_method (str): ODE の解法。 ``'Euler'`` を指定すると前進 Euler 法を用いる。
rtol (float): ``solve_ivp`` 利用時の相対誤差許容値。
atol (float): ``solve_ivp`` 利用時の絶対誤差許容値。
dense_output (bool): ``solve_ivp`` で密な出力を生成するかどうか。
matrix_type (str): 行列の種類。現在は `"symmetric"` のみ対応。
method (str): 固有分解の手法。現在は `"eigh"` のみ対応。
solver_method (str): 追跡の手法。 `'ogita_aishima'` がデフォルト。
`solve_ivp` の `method` を指定可能。
correction_method (str or None): 補正手法。 `'none'` がデフォルト。
`solver_method` で ODE ソルバー指定の場合の後処理方法を指定。
詳細は `correct_trajectory` 関数を参照。
rtol (float): `solve_ivp` 利用時の相対誤差許容値。
atol (float): `solve_ivp` 利用時の絶対誤差許容値。
dense_output (bool): `solve_ivp` で密な出力を生成するかどうか。

Returns:
EigenTrackingResults: 追跡と解析結果を含むオブジェクト。
Expand All @@ -122,16 +151,19 @@ def eigenpairtrack(

# --- Dispatch to the appropriate low-level tracker ---
try:
Qs, Lambdas, sol = _track_symmetric_eigh(
A_func,
dA_func,
t_span,
t_eval,
solver_method=solver_method,
rtol=rtol,
atol=atol,
dense_output=dense_output,
)
if solver_method == "ogita_aishima":
Qs, Lambdas, sol = _track_ogita_aishima(A_func, t_eval)
else:
Qs, Lambdas, sol = _track_symmetric_eigh(
A_func,
dA_func,
t_span,
t_eval,
solver_method=solver_method,
rtol=rtol,
atol=atol,
dense_output=dense_output,
)
success = sol.success
message = sol.message
except RuntimeError as e:
Expand Down