mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-30 16:23:41 +02:00
matrix vector product OK
This commit is contained in:
parent
3de8e31a18
commit
c5a16d0765
117
Utils/Matrix.ml
117
Utils/Matrix.ml
@ -107,7 +107,7 @@ let outer_product ?(threshold=epsilon) v1 v2 =
|
|||||||
let v1 = Vector.to_vec v1
|
let v1 = Vector.to_vec v1
|
||||||
and v2 = Vector.to_vec v2
|
and v2 = Vector.to_vec v2
|
||||||
in
|
in
|
||||||
let a = Mat.create (Vec.dim v1) (Vec.dim v2) in
|
let a = Mat.make0 (Vec.dim v1) (Vec.dim v2) in
|
||||||
ger v1 v2 a;
|
ger v1 v2 a;
|
||||||
Dense a
|
Dense a
|
||||||
| (true, false) ->
|
| (true, false) ->
|
||||||
@ -256,16 +256,71 @@ let mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b =
|
|||||||
| (Sparse a), (Sparse b) -> mmspmm transa transb a b
|
| (Sparse a), (Sparse b) -> mmspmm transa transb a b
|
||||||
|
|
||||||
|
|
||||||
|
let mv ?(sparse=false) ?(trans=`N) ?(threshold=epsilon) a b =
|
||||||
|
|
||||||
|
let f =
|
||||||
|
match trans with
|
||||||
|
| `N -> dim2
|
||||||
|
| `T -> dim1
|
||||||
|
in
|
||||||
|
if f a <> Vector.dim b then
|
||||||
|
invalid_arg "Inconsistent dimensions";
|
||||||
|
|
||||||
|
let spmv a b =
|
||||||
|
let {m ; n ; v} =
|
||||||
|
if trans = `N then
|
||||||
|
match transpose (Sparse a) with
|
||||||
|
| Sparse x -> x
|
||||||
|
| _ -> assert false
|
||||||
|
else
|
||||||
|
a
|
||||||
|
in
|
||||||
|
Array.map (fun row_a -> Vector.dot row_a b) v
|
||||||
|
|> Vec.of_array
|
||||||
|
in
|
||||||
|
|
||||||
|
let mv a b =
|
||||||
|
let f_a =
|
||||||
|
match trans with
|
||||||
|
| `N -> (fun i -> Mat.copy_row a i)
|
||||||
|
| `T -> (fun i -> Mat.col a i)
|
||||||
|
in
|
||||||
|
Vec.init (Mat.dim1 a) (fun i ->
|
||||||
|
Vector.dense_of_vec (f_a i)
|
||||||
|
|> Vector.dot b )
|
||||||
|
in
|
||||||
|
|
||||||
|
let dense_result =
|
||||||
|
match a, Vector.is_dense b with
|
||||||
|
| Dense a, true -> gemv ~trans a (Vector.to_vec b)
|
||||||
|
| Dense a, false -> mv a b
|
||||||
|
| Sparse a, true -> spmv a b
|
||||||
|
| Sparse a, false -> spmv a b
|
||||||
|
in
|
||||||
|
|
||||||
|
if sparse then
|
||||||
|
Vector.sparse_of_vec dense_result
|
||||||
|
else
|
||||||
|
Vector.dense_of_vec dense_result
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
let rec pp_matrix ppf = function
|
let rec pp_matrix ppf = function
|
||||||
| Dense m -> Util.pp_matrix ppf m
|
| Dense m -> Util.pp_matrix ppf m
|
||||||
| Sparse m -> pp_matrix ppf @@ dense_of_sparse (Sparse m)
|
| Sparse m -> pp_matrix ppf @@ dense_of_sparse (Sparse m)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
(* ---------- Unit tests ------------ *)
|
||||||
|
|
||||||
let test_case () =
|
let test_case () =
|
||||||
|
|
||||||
let d1 = 3
|
let d1 = 30
|
||||||
and d2 = 4
|
and d2 = 40
|
||||||
and d3 = 5
|
and d3 = 50
|
||||||
in
|
in
|
||||||
|
|
||||||
let x1 = Mat.map (fun x -> if abs_float x < 0.6 then 0. else x) (Mat.random d1 d2)
|
let x1 = Mat.map (fun x -> if abs_float x < 0.6 then 0. else x) (Mat.random d1 d2)
|
||||||
@ -280,6 +335,11 @@ let test_case () =
|
|||||||
and m2_s = sparse_of_mat x2
|
and m2_s = sparse_of_mat x2
|
||||||
in
|
in
|
||||||
|
|
||||||
|
let norm_diff m1 m2 =
|
||||||
|
(Mat.sub (to_mat m1) (to_mat m2)
|
||||||
|
|> Mat.syrk_trace)
|
||||||
|
in
|
||||||
|
|
||||||
let test_dimensions () =
|
let test_dimensions () =
|
||||||
Alcotest.(check int) "dim1 1" d1 (dim1 m1 );
|
Alcotest.(check int) "dim1 1" d1 (dim1 m1 );
|
||||||
Alcotest.(check int) "dim1 2" d1 (dim1 m1_s);
|
Alcotest.(check int) "dim1 2" d1 (dim1 m1_s);
|
||||||
@ -309,8 +369,8 @@ let test_case () =
|
|||||||
in
|
in
|
||||||
|
|
||||||
let test_outer () =
|
let test_outer () =
|
||||||
let x1 = Vec.init 3 (fun i -> float_of_int i)
|
let x1 = Vec.init d1 (fun i -> float_of_int i)
|
||||||
and x2 = Vec.init 4 (fun i -> float_of_int i -. 3.)
|
and x2 = Vec.init d2 (fun i -> float_of_int i -. 0.3)
|
||||||
in
|
in
|
||||||
let v1 = Vector.dense_of_vec x1
|
let v1 = Vector.dense_of_vec x1
|
||||||
and v2 = Vector.dense_of_vec x2
|
and v2 = Vector.dense_of_vec x2
|
||||||
@ -319,15 +379,43 @@ let test_case () =
|
|||||||
and v2_s = Vector.sparse_of_vec x2
|
and v2_s = Vector.sparse_of_vec x2
|
||||||
in
|
in
|
||||||
let m1 =
|
let m1 =
|
||||||
Dense (Mat.init_cols 3 4 (fun i j -> (float_of_int i) *. (float_of_int j -. 3.)))
|
Dense (Mat.init_cols d1 d2 (fun i j -> (float_of_int i) *. (float_of_int j -. 0.3)))
|
||||||
in
|
in
|
||||||
let m1_s =
|
let m1_s =
|
||||||
sparse_of_dense m1
|
sparse_of_dense m1
|
||||||
in
|
in
|
||||||
Alcotest.(check bool) "dense dense " true (m1 = outer_product v1 v2);
|
Alcotest.(check (float 1.e-10)) "dense dense " 0. (norm_diff m1 (outer_product v1 v2));
|
||||||
Alcotest.(check bool) "sparse dense " true (m1_s = outer_product v1_s v2);
|
Alcotest.(check (float 1.e-10)) "sparse dense " 0. (norm_diff m1_s (outer_product v1_s v2));
|
||||||
Alcotest.(check bool) "dense sparse" true (m1_s = outer_product v1 v2_s);
|
Alcotest.(check (float 1.e-10)) "dense sparse" 0. (norm_diff m1_s (outer_product v1 v2_s));
|
||||||
Alcotest.(check bool) "sparse sparse" true (m1_s = outer_product v1_s v2_s);
|
Alcotest.(check (float 1.e-10)) "sparse sparse" 0. (norm_diff m1_s (outer_product v1_s v2_s));
|
||||||
|
in
|
||||||
|
|
||||||
|
let test_mv () =
|
||||||
|
let y = Vec.random d2 in
|
||||||
|
let z = Vec.random d1 in
|
||||||
|
let x3 = gemv x1 y in
|
||||||
|
let x4 = gemv ~trans:`T x1 z in
|
||||||
|
let v = Vector.dense_of_vec y in
|
||||||
|
let v2 = Vector.dense_of_vec z in
|
||||||
|
let v3 = Vector.dense_of_vec x3 in
|
||||||
|
let v4 = Vector.dense_of_vec x4 in
|
||||||
|
let v_s = Vector.sparse_of_vec y in
|
||||||
|
let v2_s = Vector.sparse_of_vec z in
|
||||||
|
let norm_diff v1 v2 =
|
||||||
|
Vec.sub (Vector.to_vec v1) (Vector.to_vec v2)
|
||||||
|
|> nrm2
|
||||||
|
in
|
||||||
|
Alcotest.(check (float 1.e-10)) "dense dense 1" 0. (norm_diff (mv m1 v) v3);
|
||||||
|
Alcotest.(check (float 1.e-10)) "dense dense 2" 0. (norm_diff (mv ~trans:`T m1 v2) v4);
|
||||||
|
|
||||||
|
Alcotest.(check (float 1.e-10)) "dense sparse 3" 0. (norm_diff (mv m1 v_s) v3);
|
||||||
|
Alcotest.(check (float 1.e-10)) "dense sparse 4" 0. (norm_diff (mv ~trans:`T m1 v2_s) v4);
|
||||||
|
|
||||||
|
Alcotest.(check (float 1.e-10)) "sparse dense 5" 0. (norm_diff (mv m1_s v) v3);
|
||||||
|
Alcotest.(check (float 1.e-10)) "sparse dense 6" 0. (norm_diff (mv ~trans:`T m1_s v2) v4);
|
||||||
|
|
||||||
|
Alcotest.(check (float 1.e-10)) "sparse sparse 7" 0. (norm_diff (mv m1_s v_s) v3);
|
||||||
|
Alcotest.(check (float 1.e-10)) "sparse sparse 8" 0. (norm_diff (mv ~trans:`T m1_s v2_s) v4);
|
||||||
in
|
in
|
||||||
|
|
||||||
let test_mm () =
|
let test_mm () =
|
||||||
@ -339,10 +427,6 @@ let test_case () =
|
|||||||
and m5 = dense_of_mat x2 |> transpose
|
and m5 = dense_of_mat x2 |> transpose
|
||||||
and m5_s = sparse_of_mat x2 |> transpose
|
and m5_s = sparse_of_mat x2 |> transpose
|
||||||
in
|
in
|
||||||
let norm_diff m1 m2 =
|
|
||||||
(Mat.sub (to_mat m1) (to_mat m2)
|
|
||||||
|> Mat.syrk_trace)
|
|
||||||
in
|
|
||||||
Alcotest.(check (float 1.e-10)) "dense dense 1" 0. (norm_diff (mm m1 m2) m3);
|
Alcotest.(check (float 1.e-10)) "dense dense 1" 0. (norm_diff (mm m1 m2) m3);
|
||||||
Alcotest.(check (float 1.e-10)) "dense dense 2" 0. (norm_diff (mm ~transa:`T m4 m2) m3);
|
Alcotest.(check (float 1.e-10)) "dense dense 2" 0. (norm_diff (mm ~transa:`T m4 m2) m3);
|
||||||
Alcotest.(check (float 1.e-10)) "dense dense 3" 0. (norm_diff (mm ~transb:`T m1 m5) m3);
|
Alcotest.(check (float 1.e-10)) "dense dense 3" 0. (norm_diff (mm ~transb:`T m1 m5) m3);
|
||||||
@ -367,7 +451,8 @@ let test_case () =
|
|||||||
"Conversion", `Quick, test_conversion;
|
"Conversion", `Quick, test_conversion;
|
||||||
"Dimensions", `Quick, test_dimensions;
|
"Dimensions", `Quick, test_dimensions;
|
||||||
"Transposition", `Quick, test_transpose;
|
"Transposition", `Quick, test_transpose;
|
||||||
"Oueter product", `Quick, test_outer;
|
"Outer product", `Quick, test_outer;
|
||||||
|
"Matrix Vector", `Quick, test_mv;
|
||||||
"Matrix Matrix", `Quick, test_mm;
|
"Matrix Matrix", `Quick, test_mm;
|
||||||
]
|
]
|
||||||
|
|
||||||
|
71
Utils/Matrix.mli
Normal file
71
Utils/Matrix.mli
Normal file
@ -0,0 +1,71 @@
|
|||||||
|
(* Sparse or dense matrices *)
|
||||||
|
|
||||||
|
open Lacaml.D
|
||||||
|
|
||||||
|
type t
|
||||||
|
|
||||||
|
(** {1 Accessors} *)
|
||||||
|
|
||||||
|
val is_sparse : t -> bool
|
||||||
|
(** True is the matrix is sparse. *)
|
||||||
|
|
||||||
|
val is_dense : t -> bool
|
||||||
|
(** True is the matrix is dense. *)
|
||||||
|
|
||||||
|
val dim1 : t -> int
|
||||||
|
(** First dimension of the matrix *)
|
||||||
|
|
||||||
|
val dim2 : t -> int
|
||||||
|
(** Secons dimension of the matrix *)
|
||||||
|
|
||||||
|
val get : t -> int -> int -> float
|
||||||
|
(** [get m i j] returns {% $\mathbf{m}_{ij}$ %}. *)
|
||||||
|
|
||||||
|
|
||||||
|
(** {1 Converters } *)
|
||||||
|
|
||||||
|
val to_mat : t -> Mat.t
|
||||||
|
(** Convert into a Lacaml Mat. *)
|
||||||
|
|
||||||
|
val to_vector_array : t -> Vector.t array
|
||||||
|
(** Convert the matrix into an array of column vectors. *)
|
||||||
|
|
||||||
|
val sparse_of_dense : ?threshold:float -> t -> t
|
||||||
|
(** Creates a sparse matrix from a dense matrix. Default threshold is {!Constants.epsilon}. *)
|
||||||
|
|
||||||
|
val dense_of_sparse : t -> t
|
||||||
|
(** Creates a dense matrix from a sparse matrix. *)
|
||||||
|
|
||||||
|
val dense_of_mat : Mat.t -> t
|
||||||
|
(** Create a dense matrix from a Lacaml Mat *)
|
||||||
|
|
||||||
|
val sparse_of_mat : ?threshold:float -> Lacaml.D.Mat.t -> t
|
||||||
|
(** Create a sparse matrix from a Lacaml Mat. Default threshold is {!Constants.epsilon}. *)
|
||||||
|
|
||||||
|
val sparse_of_vector_array : Vector.t array -> t
|
||||||
|
(** Create a sparse matrix from an array of column vectors. *)
|
||||||
|
|
||||||
|
val transpose : t -> t
|
||||||
|
(** Returns the transposed matrix. *)
|
||||||
|
|
||||||
|
|
||||||
|
(** {1 Operations} *)
|
||||||
|
|
||||||
|
val outer_product : ?threshold:float -> Vector.t -> Vector.t -> t
|
||||||
|
(** Creates the matrix formed by the outer product of two vectors. *)
|
||||||
|
|
||||||
|
val mm : ?transa:trans3 -> ?transb:trans3 -> ?threshold:float -> t -> t -> t
|
||||||
|
(** Matrix multiplication *)
|
||||||
|
|
||||||
|
val mv : ?sparse:bool -> ?trans:trans3 -> ?threshold:float -> t -> Vector.t -> Vector.t
|
||||||
|
(** Matrix Vector product *)
|
||||||
|
|
||||||
|
|
||||||
|
(** {1 Printers } *)
|
||||||
|
|
||||||
|
val pp_matrix : Format.formatter -> t -> unit
|
||||||
|
|
||||||
|
|
||||||
|
(** {1 Unit testing} *)
|
||||||
|
|
||||||
|
val test_case : unit -> (string * [> `Quick ] * (unit -> unit)) list
|
@ -1,7 +1,7 @@
|
|||||||
open Lacaml.D
|
|
||||||
|
|
||||||
(* Sparse or dense vectors *)
|
(* Sparse or dense vectors *)
|
||||||
|
|
||||||
|
open Lacaml.D
|
||||||
|
|
||||||
type t
|
type t
|
||||||
|
|
||||||
(** {1 Accessors} *)
|
(** {1 Accessors} *)
|
||||||
|
Loading…
Reference in New Issue
Block a user