defmodule Chi2fit.Matrix do
# Copyright 2016-2017 Pieter Rijken
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
@moduledoc """
This module provides matrix inverse operations and supporting functions.
It provides 2 types of matrix norms and an iterative approach to calculating the matrix inverse.
The implementation is based on the work [1].
## References
[1] F. Soleymani, A Rapid Numerical Algorithm to Compute Matrix Inversion, International Journal of Mathematics and
Mathematical Sciences, Volume 2012, Article ID 134653, doi:10.1155/2012/134653
"""
@inverse_tolerance 1.0e-8
@default_inverse_iterations 500
@default_range 100
@default_size 100
@default_algorithm :lie
@typedoc """
A list of numbers
"""
@type vector :: [number]
@typedoc """
A list of vectors (list of lists of numbers)
"""
@type matrix :: [vector]
@doc """
Constructs a unit matrix of size n. All diagonal elements have value 1 and the rest has value 0.
## Examples
iex> unit(3)
[[1, 0, 0], [0, 1, 0], [0, 0, 1]]
iex> unit(0)
** (ArgumentError) Illegal argument '0'
iex> unit -1
** (ArgumentError) Illegal argument '-1'
iex> unit 1.3
** (ArgumentError) Illegal argument '1.3'
"""
@spec unit(n :: pos_integer) :: [[0|1]]
def unit(n) when is_integer(n) and n > 0 do
{result,_} = List.duplicate(0,n) |> List.duplicate(n) |> Enum.reduce({[],-1}, fn (list,{acc,m}) -> {[list |> List.replace_at(m,1)|acc],m-1} end)
result
end
def unit(n), do: raise ArgumentError, message: "Illegal argument '#{inspect n}'"
defp abssum(list), do: list |> Enum.map(&abs/1) |> Enum.sum
@doc """
Calculates the norm of the matrix as the sum of the absolutes value of all elements.
## Example
iex> norm [[1,2,3],[4,5,6],[7,8,9]]
45
"""
@spec norm(matrix) :: number
def norm(matrix), do: matrix |> Enum.map(&abssum/1) |> Enum.sum
@doc """
Calculates the norm of the matrix. All absolute values of the elements of each row are summed. The maximum
value is returned
## Example
iex> norm_1 [[1,2,3],[4,5,6],[7,8,9]]
24
"""
@spec norm_1(matrix) :: number
def norm_1(matrix), do: matrix |> Enum.map(&abssum/1) |> Enum.max
@doc """
Calculates the norm of the matrix as `norm_1/1` but over the columns instead of over the rows.
## Example
iex> norm_inf [[1,2,3],[4,5,6],[7,8,9]]
18
"""
@spec norm_inf(matrix) :: number
def norm_inf(matrix), do: matrix |> transpose |> norm_1
defmacrop telescope(matrix, []), do: quote(do: unit(length unquote(matrix)))
defmacrop telescope(matrix, [a|rest]) do
quote do
mat = unquote(matrix)
subtract(scalar_multiply(unit(length mat),unquote(a)),multiply(mat,telescope(mat, unquote(rest))))
end
end
defp findv0(:way2, matrix, _options) do
v0 = matrix |> transpose |> scalar_multiply(1.0/norm_1(matrix)/norm_inf(matrix))
test = matrix |> length |> unit |> subtract(multiply(matrix,v0)) |> norm_1
{v0,test}
end
defp findv0(:way3, matrix, options) do
range = options[:range] || @default_range
size = options[:size] || @default_size
List.duplicate(0,size)
|> Enum.map(fn (_x)->range*(2*:rand.uniform() - 1) end)
|> Enum.reduce({nil,:infinity},fn
(factor,{_,:infinity}) ->
v0 = matrix |> length |> unit |> scalar_multiply(factor)
test = matrix |> length |> unit |> subtract(multiply(matrix,v0)) |> norm_1 |> abs
{v0,test}
(factor,{v,error}) ->
v0 = matrix |> length |> unit |> scalar_multiply(factor)
test = matrix |> length |> unit |> subtract(multiply(matrix,v0)) |> norm_1 |> abs
if test < error, do: {v0,test}, else: {v,error}
end)
end
@doc """
Returns the matrix inverse of the argument.
## Options
`:tolerance` - Iterate until the `norm_1/1` of I-AV is less than this value
`:algorithm` - Four algorithms are supported: `:hotelling_bodewig` (second order), `:lie` (third order),
`:krishnamurthy_sen` (sixth order), and `:soleymani` (seventh order); defaults to `#{inspect @default_algorithm}`
`:max_iterations` - Maximum number of iterations to perform; defaults to #{@default_inverse_iterations}
`:range` - Range of values from -range to +range as a multiple of the unit matrix to try as an estimate
of the inverse matric; defaults to #{@default_range}
`:size` - Number of tries to estimate initial inverse; defautls to #{@default_size}
## Examples
iex> inverse [[3]]
{:ok,[[0.3333333333333333]]}
iex> inverse [[1,2],[3,4]]
{:ok,[[-2.0, 1.0], [1.5, -0.5]]}
iex> inverse([[3,2,0],[0,0,1],[2,-2,1]]) |> elem(1) |> Enum.map(fn row -> Enum.map(row, & Float.round(&1,10)) end)
[[0.2, -0.2, 0.2], [0.2, 0.3, -0.3], [0.0, 1.0, 0.0]]
iex> inverse([[3,2,0],[0,0,1],[2,-2,1]], algorithm: :soleymani) |> elem(1) |> Enum.map(fn row -> Enum.map(row, & Float.round(&1,14)) end)
[[0.2, -0.2, 0.2], [0.2, 0.3, -0.3], [0.0, 1.0, 0.0]]
iex> inverse([[3,2,0],[0,0,1],[2,-2,1]], tolerance: 1.0e-15) |> elem(1) |> Enum.map(fn row -> Enum.map(row, & Float.round(&1,14)) end)
[[0.2, -0.2, 0.2], [0.2, 0.3, -0.3], [0.0, 1.0, 0.0]]
For matrices that have no inverse:
iex> try do inverse [[1,2,3],[4,5,6],[7,8,9]] catch x->x end
:no_inverse
"""
@spec inverse(matrix, options :: Keyword.t) :: {:ok,inverse :: matrix} | :failed_to_find_v0 | :no_inverse | {:failed_to_reach_tolerance,inverse :: matrix,error :: float}
def inverse(matrix, options \\ [])
def inverse([[x]], _options), do: {:ok,[[1.0/x]]}
def inverse([[x1,x2],[y1,y2]], _options) when x1*y2-x2*y1==0.0, do: throw {:impossible_inverse,"Zero determinant"}
def inverse([[x1,x2],[y1,y2]], _options), do: {:ok,[[y2,-x2],[-y1,x1]] |> scalar_multiply(1.0/(x1*y2-x2*y1))}
def inverse(matrix, options) do
max_iter = options[:max_iterations] || @default_inverse_iterations
{v0,test} = findv0(:way2,matrix, options)
try do
if test < 2.0 do
iterate(matrix,v0,max_iter,options)
else
{v0,test} = findv0(:way3,matrix,options)
if test < 1.0 do
iterate(matrix,v0,max_iter,options)
else
throw :no_v0
end
end
rescue
_e in ArithmeticError ->
throw :no_inverse
catch
{:impossible_inverse,v,error} ->
throw {:failed_to_reach_tolerance,v,error}
:no_v0 ->
throw :failed_to_find_v0
end
end
defp forward(:hotelling_bodewig, av), do: telescope(av,[2.0])
defp forward(:lie, av), do: telescope(av,[3.0,3.0])
defp forward(:krishnamurthy_sen, av), do: telescope(av,[2.0]) |> multiply(telescope(av,[3.0,3.0]) |> multiply(telescope(av,[1.0,1.0])))
defp forward(:soleymani, av), do: telescope(av,[120.0,393.0,735.0,861.0,651.0,315.0,93.0,15.0]) |> scalar_multiply(1/16.0)
defp iterate(matrix,v0,0,_options) do
throw {:impossible_inverse,v0,subtract(unit(length matrix),multiply(matrix,v0)) |> norm_1}
end
defp iterate(matrix,v0,max,options) when is_integer(max) and max > 0 do
tolerance = options[:tolerance] || @inverse_tolerance
algorithm = options[:algorithm] || @default_algorithm
u = unit(length matrix)
av = multiply(matrix,v0)
test = subtract(u,av) |> norm_1
unless test < tolerance do
matrix |> iterate(multiply(v0,forward(algorithm, av)),max-1,options)
else
{:ok,v0}
end
end
@doc """
Returns the diagonal elements of the matrix as a vector.
## Example
iex> diagonal [[1,2,3],[4,5,6],[7,8,9]]
[1, 5, 9]
"""
@spec diagonal(matrix) :: vector
def diagonal(matrix) do
matrix |> Enum.reduce({[],0}, fn (row, {acc,index})->{[Enum.at(row,index)|acc],index+1} end) |> elem(0) |> Enum.reverse
end
@doc """
Returns a matrix with the supplied vector as its diagonal elements.
## Examples
iex> from_diagonal [1,5,9]
[[1, 0, 0], [0, 5, 0], [0, 0, 9]]
"""
@spec from_diagonal(vector) :: matrix
def from_diagonal(vector) do
vector |> Enum.reduce({[],0}, fn (elem, {acc,index})->{[List.duplicate(0,length(vector)) |> List.replace_at(index,elem)|acc],index+1} end) |> elem(0) |> Enum.reverse
end
@doc """
Calculates the inner product of two vectors.
## Examples
iex> dotproduct [1,2], [3,4]
11
iex> dotproduct [], []
0
iex> dotproduct [1,2], []
** (ArgumentError) Vectors of unequal length
iex> dotproduct [1,2], [1]
** (ArgumentError) Vectors of unequal length
"""
@spec dotproduct(vector,vector) :: number
def dotproduct(vector1, vector2, sum \\ 0)
def dotproduct(vector1, vector2, _sum) when length(vector1) != length(vector2) do
raise ArgumentError, message: "Vectors of unequal length"
end
def dotproduct([], [], sum), do: sum
def dotproduct([v1|rest1], [v2|rest2], sum), do: dotproduct(rest1, rest2, sum + v1*v2)
@doc """
Subtracts two matrices and returns the result.
"""
@spec subtract(matrix,matrix) :: matrix
def subtract(matrix1, matrix2), do: subtract(matrix1,matrix2,[])
defp subtract([], [], result), do: Enum.reverse(result)
defp subtract([row1|rest1], [row2|rest2], result) do
subtract(rest1, rest2, [subtractv(row1,row2)|result])
end
defp subtractv(vector1, vector2, result \\ [])
defp subtractv([], [], result), do: Enum.reverse(result)
defp subtractv([v1|rest1], [v2|rest2], result) do
subtractv(rest1, rest2, [v1-v2|result])
end
defp scalar_multiply(matrix, scalar) when is_number(scalar) do
matrix |> Enum.map(fn row -> Enum.map(row, & &1*scalar) end)
end
defp multiply(matrix1, matrix2, result \\ [])
defp multiply([], _matrix2, result), do: Enum.reverse(result)
defp multiply([row|rest1], matrix2, result) do
multiply(rest1, matrix2, [matrix2
|> Stream.zip
|> Enum.reduce([], fn col, acc -> [dotproduct(row, Tuple.to_list(col))|acc] end)
|> Enum.reverse | result])
end
@doc """
Returns the tranpose of the matrix
## Examples:
iex> transpose [ [1] ]
[[1]]
iex> transpose [ [1,2], [3,4] ]
[[1, 3], [2, 4]]
"""
@spec transpose(matrix) :: matrix
def transpose(matrix), do: matrix |> Stream.zip |> Stream.map(&Tuple.to_list/1) |> Enum.into([])
@doc """
Adds two vectors.
## Examples
iex> add [1,2], [3,4]
[4,6]
iex> add [], []
[]
iex> add [1], [5]
[6]
"""
@spec add(vector, vector) :: vector
def add(vector1, vector2), do: Stream.zip(vector1, vector2) |> Enum.map(fn {v1,v2} -> v1 + v2 end)
@doc """
Calculates the determinant of the matrix.
"""
@spec det(matrix) :: number
def det([x]), do: x
def det([ [a11,a12], [a21,a22]] ), do: a11*a22 - a12*a21
end
