10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-07 14:43:41 +01:00
QCaml/Utils/DIIS.mli

63 lines
1.9 KiB
OCaml
Raw Normal View History

2018-05-31 10:47:31 +02:00
(** Direct Inversion of the Iterative Subspace algorithm.
2018-05-31 10:56:01 +02:00
At iteration {% $m$ %}, one has:
2018-05-31 10:47:31 +02:00
2018-05-31 10:56:01 +02:00
- {% $\mathbf{p}_m$ %}, a vector of parameters
- {% $\mathbf{e}_m$ %}, an approximate error vector
The DIIS approximate solution for iteration $m+1$ is given by
2018-05-31 10:47:31 +02:00
{% \begin{align*}
2018-05-31 10:56:01 +02:00
\mathbf{p}_{m+1} & = \sum_{i=1}^m c_i (\mathbf{p}^f + \mathbf{e}_i) \\
& = \sum_{i=1}^m c_i \mathbf{p}^f + \sum_i c_i \mathbf{e}_i
2018-05-31 10:47:31 +02:00
\end{align*} %}
2018-05-31 10:56:01 +02:00
where {% $\mathbf{p}^f$ %} is the exact solution. One wants to minimize the
norm of the error vector imposing the constraint that {% $\sum_{i=1}^m c_i = 1$ %}
with a Langrange multiplier {% $\lambda$ %}.
2018-05-31 10:47:31 +02:00
{%
\begin{align*}
2018-05-31 10:56:01 +02:00
\mathcal{L} & = ||\sum_i c_i \mathbf{e}_i||^2 - \lambda \left(\sum_i c_i - 1\right) \\
2018-06-27 13:13:59 +02:00
& = \sum_{ij} c_i c_j B_{ij} - \lambda \left(\sum_{i=1}^m c_i - 1\right)
2018-05-31 10:47:31 +02:00
\end{align*}
2018-05-31 10:56:01 +02:00
%}
with
{% $B_{ij} = \langle \mathbf{e}_i | \mathbf{e}_j \rangle$ %}.
2018-05-31 10:47:31 +02:00
Equating zero to the derivatives of {% $\mathcal{L}$ %} with respect to {% $c_i$ %} and {% $\lambda$ %} leads to
{% \begin{equation*}
\begin{bmatrix}
2018-05-31 16:46:45 +02:00
B_{11} & B_{12} & B_{13} & ... & B_{1m} & 1 \\
B_{21} & B_{22} & B_{23} & ... & B_{2m} & 1 \\
B_{31} & B_{32} & B_{33} & ... & B_{3m} & 1 \\
2018-05-31 10:47:31 +02:00
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
2018-05-31 16:46:45 +02:00
B_{m1} & B_{m2} & B_{m3} & ... & B_{mm} & 1 \\
2018-05-31 10:47:31 +02:00
1 & 1 & 1 & ... & 1 & 0
2018-05-31 16:46:45 +02:00
\end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_m \\ -\lambda \end{bmatrix}=
2018-05-31 10:47:31 +02:00
\begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}
\end{equation*}
%}
The coefficients are then used to update {% $\mathbf{p}$ %} as
{% $$
\mathbf{p}_{m+1}=\sum_{i=1}^m c_i\mathbf{p}_i.
$$ %}
*)
type t
2018-05-31 16:46:45 +02:00
val make : ?mmax:int -> unit -> t
(** Initialize DIIS with a maximum size.*)
2018-05-31 10:47:31 +02:00
2018-05-31 16:46:45 +02:00
val append : p:Lacaml.D.Vec.t -> e:Lacaml.D.Vec.t -> t -> t
2018-05-31 10:47:31 +02:00
(** Append a parameter vector [p] and the corresponding error vector [e]. *)
val next : t -> Lacaml.D.Vec.t
(** Returns a new parameter vector. *)