逆行列の計算

必要になって何となく OCaml で書いて何となく晒す。

なんとなく面倒臭そうなイメージを持ってたんですが，実際に書いてみたらそうでもありませんでした。

破壊的代入の塊であるという点では関数型らしくありませんが，配列の破壊的操作が Array.iteri で思いの外すっきり書けたあたりは関数型言語の特性かなあ，と思いました。

let init_matrix m n f =
  Array.init m (fun i -> Array.init n (fun j -> f i j))

let make_unit_matrix n =
  init_matrix n n (fun i j -> if i = j then 1.0 else 0.0)

(* destructive version of Array.map; works as (map-into a f a) in CL *)
let amap f a =
  Array.iteri (fun i x -> a.(i) <- f x) a

(* another variant: (map-into a f a b) *)
let amap2 f a b =
  Array.iteri (fun i x -> a.(i) <- f x b.(i)) a

(* some destructive operations on rows *)
let swap_rows m i j =
  let ri, rj = m.(i), m.(j) in
    m.(i) <- rj; m.(j) <- ri

let subtract_row m i j x =
  amap2 (fun a b -> a -. b *. x) m.(i) m.(j)

let magnify_row m i x =
  amap (( *. ) x) m.(i)

(* sweeping algorithm *)
let rec find_pivot m i' i =
  if m.(i').(i) <> 0.0 then i' else find_pivot m (i' + 1) i

let sweep_ m i =
  (let i' = find_pivot m i i in swap_rows m i i');
  magnify_row m i (1.0 /. m.(i).(i));
  for i' = 0 to Array.length m - 1 do
    if i <> i' then subtract_row m i' i (m.(i').(i))
  done

let sweep m = for i = 0 to Array.length m - 1 do sweep_ m i done; m

(* ad-hoc stuff *)
let concatenate_matrix m1 m2 =
  Array.init (Array.length m1) (fun i -> Array.append (m1.(i)) (m2.(i)))

let extract_rhs m =
  let n = Array.length m in
    Array.init n (fun i -> Array.sub (m.(i)) n n)

(* main function *)
let inverse_matrix m =
  let e = make_unit_matrix (Array.length m) in
  let m' = concatenate_matrix m e in
    extract_rhs (sweep m')