Improved sparse vectors

Utils/Davidson.ml

@@ -6,7 +6,7 @@ let make
     ?guess
     ?(n_states=8)
     ?(n_iter=10)
-    ?(threshold=1.e-8)
+    ?(threshold=1.e-6)
     diagonal
     matrix_prod
   =
@@ -24,16 +24,16 @@ let make
   let random_vectors =
     let random_vector k =
       Vec.init n (fun i ->
-          let r1 = Random.float 1.
-          and r2 = Random.float 1.
-          in
-          let a = sqrt (-2. *. log r1)
-          and b = Constants.two_pi *. r2
-          in
           if i<k then 0.
-          else if i>k then
+          else if i=k then 1.e5
-            a *. cos b
+          else
-          else 100.0
+            let r1 = Random.float 1.
+            and r2 = Random.float 1.
+            in
+            let a = sqrt (-2. *. log r1)
+            and b = Constants.two_pi *. r2 in
+            let c = a *. cos b in
+            if abs_float c > 1.e-1 then c else 0.
         )
       |> Util.normalize
     in
@@ -73,7 +73,7 @@ let make
       matrix_prod (
         u_new_ortho
         |> Mat.of_col_vecs_list
-        |> Matrix.dense_of_mat )
+        |> Matrix.sparse_of_mat )
       |> Matrix.to_mat
       |> Mat.to_col_vecs_list
     in             

Utils/Vector.ml

@@ -2,10 +2,16 @@ open Lacaml.D
 
 let epsilon = Constants.epsilon
 
+type index_value = 
+  {
+    index: int;
+    value: float
+  }
+
 type sparse_vector =
   {
     n: int;
-    v: (int*float) list
+    v: index_value list
   }
 
 type t =
@@ -23,13 +29,22 @@ let is_dense = function
   | Dense  _ -> true
 
 
+exception Found of float
+
 let get = function
   | Dense  v ->  (fun i -> v.{i})
   | Sparse { n ; v } ->  (fun i ->
       if i < 1 || i > n then invalid_arg "index out of bounds";
-      match List.assoc_opt i v with
+      try
-      | Some x -> x
+        List.iter (fun {index ; value} ->
-      | None   -> 0. )
+            if index=i then
+              raise (Found value)) v; 
+        raise Not_found
+      with
+      | Not_found -> 0.
+      | Found x -> x
+    )
+
 
 
 let dim = function
@@ -47,14 +62,14 @@ let sparse_of_dense ?(threshold=epsilon) = function
         if abs_float x < threshold then
           aux accu (i-1)
         else
-          aux ((i, x)::accu) (i-1)
+          aux ({index=i ; value=x}::accu) (i-1)
     in
     let n = Vec.dim v in
     Sparse { n ; v=aux [] n }
 
 
 let rec to_assoc_list ?(threshold=epsilon) = function
-  | Sparse {n ; v} -> v
+  | Sparse {n ; v} -> List.map (fun {index ; value} -> (index, value)) v
   | Dense v -> to_assoc_list @@ sparse_of_dense ~threshold (Dense v)
 
 
@@ -62,7 +77,7 @@ let dense_of_sparse = function
   | Dense _ -> invalid_arg "Expected a sparse vector"
   | Sparse {n ; v} -> 
     let v' = Vec.make0 n in
-    List.iter (fun (i, x) -> v'.{i} <- x) v;
+    List.iter (fun {index ; value} -> v'.{index} <- value) v;
     Dense v'
 
 
@@ -74,7 +89,10 @@ let sparse_of_vec ?(threshold=epsilon) v =
   |> sparse_of_dense ~threshold 
 
 
-let sparse_of_assoc_list n v = Sparse { n ; v }
+let sparse_of_assoc_list n v =
+  Sparse { n ;
+    v = List.map (fun (index, value) -> {index ; value}) v
+  }
 
 
 let rec to_vec = function
@@ -85,67 +103,22 @@ let rec to_vec = function
 
 let scale ?(threshold=epsilon) x = function
   | Dense  v -> let v' = copy v in (scal x v'; Dense v')
-  | Sparse {n ; v} -> Sparse {n ; v=List.map (fun (i,y) -> let z = x *. y in
+  | Sparse {n ; v} ->
-     if abs_float z > threshold then Some (i, z) else None ) v |> Util.list_some }
+    Sparse {n ; v=List.map (fun {index ; value} ->
+        let z = x *. value in
+        if abs_float z > threshold then
+          Some {index ; value=z}
+        else
+          None
+      ) v |> Util.list_some }
 
 
 let rec neg = function
   | Dense  v -> Dense (Vec.neg v)
-  | Sparse {n ; v} -> Sparse {n ; v=List.map (fun (i,y) -> (i, -. y)) v}
+  | Sparse {n ; v} ->
+      Sparse {n ; v=List.map (fun {index ; value} -> {index ; value = -. value}) v}
 
 
-(*
-let rec add ?(threshold=epsilon) x y = 
-  if dim x <> dim y then
-    invalid_arg "Inconsistent dimensions";
-
-  match x, y with
-  | Dense  x      , Dense  y -> Dense (Vec.add x y)
-  | Sparse {n ; v}, Dense  y ->
-    let v' = copy y in
-    List.iter (fun (i, x) -> v'.{i} <- v'.{i} +. x) v;
-    sparse_of_vec ~threshold v'
-  | Sparse {n ; v}, Sparse {n=n' ; v=v'} -> 
-    begin
-      let rec aux accu v1 v2 =
-        match v1, v2 with
-        | [], [] -> {n ; v=List.rev accu}
-        | ((i, x)::v1), []  -> 
-          aux ((i, x)::accu) v1 []
-        | [], ((j, y)::v2) -> 
-          aux ((j, y)::accu) [] v2
-        | ((i, x)::v1), ((j, y)::v2) -> 
-          if i = j then
-            begin
-              let z = x +. y in
-              if abs_float z > threshold then
-                aux ((i, (x +. y))::accu) v1 v2
-              else
-                aux accu v1 v2
-            end
-          else if i < j then
-            begin
-              if abs_float x > threshold then
-                aux ((i, x)::accu) v1 ((j, y)::v2)
-              else
-                aux accu v1 ((j, y)::v2)
-            end
-          else
-            begin
-              if abs_float y > threshold then
-                aux ((j, y)::accu) ((i, x)::v1) v2
-              else
-                aux accu ((i, x)::v1) v2
-            end
-      in
-      Sparse (aux [] v v')
-    end
-  | x, y -> add ~threshold y x
-
-
-let sub ?(threshold=epsilon) x y = add ~threshold x (neg y)
-*)
-
 
 let axpy ?(threshold=epsilon) ?(alpha=1.) x y = 
   if dim x <> dim y then
@@ -156,47 +129,51 @@ let axpy ?(threshold=epsilon) ?(alpha=1.) x y =
   | Sparse {n ; v}, Dense  y ->
     begin
       let v' = copy y in
-      List.iter (fun (i, x) -> v'.{i} <- v'.{i} +. alpha *. x) v;
+      List.iter (fun {index ; value} -> v'.{index} <- v'.{index} +. alpha *. value) v;
       sparse_of_vec ~threshold v'
     end
   | Dense  x      , Sparse {n ; v}  ->
     begin
       let v' = copy x in
       scal alpha v';
-      List.iter (fun (i, y) -> v'.{i} <- v'.{i} +. y) v;
+      List.iter (fun {index ; value} -> v'.{index} <- v'.{index} +. value) v;
       sparse_of_vec ~threshold v'
     end
   | Sparse {n ; v}, Sparse {n=n' ; v=v'} -> 
     begin
       let rec aux accu v1 v2 =
         match v1, v2 with
-        | []          , []           -> {n ; v=List.rev accu}
+        | ({index=i ; value=x}::r1), ({index=j ; value=y}::r2) -> 
-        | ((i, x)::v1), []           -> aux ((i, x)::accu) v1 []
+        begin
-        | []          , ((j, y)::v2) -> aux ((j, y)::accu) [] v2
+          match compare i j with
-        | ((i, x)::v1), ((j, y)::v2) -> 
+          | -1 -> 
-          if i = j then
-            begin
-              let z = alpha *. x +. y in
-              if abs_float z > threshold then
-                aux ((i, z)::accu) v1 v2
-              else
-                aux accu v1 v2
-            end
-          else if i < j then
             let z = alpha *. x in
-            begin
+            let new_accu = 
               if abs_float z > threshold then
-                aux ((i, z)::accu) v1 ((j, y)::v2)
+                {index=i ; value=z} :: accu
               else
-                aux accu v1 ((j, y)::v2)
+                accu
-            end
+            in aux new_accu r1 v2
-          else
+          | 1 ->
-            begin
+            let new_accu = 
               if abs_float y > threshold then
-                aux ((j, y)::accu) ((i, x)::v1) v2
+                {index=j ; value=y} :: accu
               else
-                aux accu ((i, x)::v1) v2
+                accu
-            end
+            in aux new_accu v1 r2
+          | 0 ->
+            let z = alpha *. x +. y in
+            let new_accu =
+              if abs_float z > threshold then
+                {index=i ; value=z} :: accu
+              else
+                accu
+            in aux new_accu r1 r2
+          | _ -> assert false
+          end
+        | ({index=i ; value=x}::r1), []  -> aux ({index=i ; value=x}::accu) r1 []
+        | [] , ({index=j ; value=y}::r2) -> aux ({index=j ; value=y}::accu) [] r2
+        | [] , [] -> {n ; v=List.rev accu}
       in
       Sparse (aux [] v v')
     end
@@ -210,7 +187,7 @@ let pp_vector ppf = function
   | Sparse {n ; v} -> 
     begin
       Format.fprintf ppf "@[[ %d | " n;
-      List.iter (fun (i,x) -> Format.fprintf ppf "@[(%d, %f); @]" i x) v;
+      List.iter (fun {index ; value} -> Format.fprintf ppf "@[(%d, %f); @]" index value) v;
       Format.fprintf ppf "]@]"
     end
   
@@ -225,18 +202,26 @@ let dot v v' =
   let d_sp v' {n ; v} =
     if n <> Vec.dim v' then
       invalid_arg "Inconsistent dimensions";
-    List.fold_left (fun accu (i, v_i) -> accu +. v_i *. v'.{i}) 0. v
+    List.fold_left (fun accu {index ; value} -> accu +. value *. v'.{index}) 0. v
   in
 
   let sp_sp {n ; v} {n=n' ; v=v'} =
     if n <> n' then
       invalid_arg "Inconsistent dimensions";
 
-    List.fold_left (fun accu (i, v_i) -> 
+    let rec aux accu = function
-        match List.assoc_opt i v' with
+      | (({index=i ; value=v1} :: r1) as s1), (({index=j ; value=v2}::r2) as s2)-> 
-        | Some w_i -> accu +. v_i *. w_i
+        begin
-        | None     -> accu 
+          match compare i j with
-      ) 0. v
+          | -1 -> aux  accu  (r1, s2)
+          |  1 -> aux  accu  (s1, r2)
+          |  0 -> aux (accu +. v1 *. v2) (r1, r2)
+          |  _ -> assert false
+        end
+      | ([], _ ) 
+      | (_ , []) -> accu
+    in
+    aux 0. (v, v')
   in
 
   match v, v' with 

Utils/Vector.mli

@@ -60,7 +60,7 @@ val sub : ?threshold:float -> t -> t -> t
 (** Subtract two vectors *)
 
 val axpy : ?threshold:float -> ?alpha:float -> t -> t -> t
-(** $a \mathfb{x} + \mathfb{y}$ *)
+(** {% $a \mathbf{x} + \mathbf{y}$ %} *)
 
 val dot : t -> t -> float
 (** Dot product. *)