Current section

25 Versions

Jump to

Compare versions

4 files changed
+343 additions
-303 deletions
  @@ -1,6 +1,6 @@
1 1 # elixir-matrix-operation
2 2 Matrix operation library in Elixir.
3 - It is described the brief explanation (numerical_formula.pdf) for a mathematical description.
3 + It is described the brief explanation (documents/latex_out/numerical_formula.pdf) for a mathematical description.
4 4
5 5 ## Notice
6 6 The column and row numbers are specified by a integer from 1 (not from 0).
  @@ -16,7 +16,7 @@ You can install this package by adding this code to dependencies in your mix.exs
16 16 ```elixir
17 17 def deps do
18 18 [
19 - {:matrix_operation, "~> 0.3.4"}
19 + {:matrix_operation, "~> 0.3.5"}
20 20 ]
21 21 end
22 22 ```
  @@ -119,22 +119,15 @@ MatrixOperation.tensor_product([[3, 2, 3], [2, 1, 2]], [[2, 3, 1], [2, 1, 2], [3
119 119 ]
120 120 ]
121 121 ```
122 - * Eigenvalue (Algebraic method for 2×2 or 3×3 matrix)
122 + * Calculate eigenvalue by algebra method (Algebraic method for 2×2 or 3×3 matrix)
123 123 ```elixir
124 - MatrixOperation.eigenvalue([[6, -3], [4, -1]])
124 + MatrixOperation.eigenvalue_algebra([[6, -3], [4, -1]])
125 125 [3.0, 2.0]
126 126 ```
127 - * Singular value (2×n or n×2 or 3×n or n×3 matrix)
127 + * Singular value by QR decomposition
128 128 ```elixir
129 - MatrixOperation.singular_value([[0, 1], [1, 0], [1, 0]])
130 - [1.4142135623730951, 1.0]
131 - MatrixOperation.singular_value([[2, 2, 2, 2], [1, -1, 1, -1], [-1, 1, -1, 1]])
132 - [4.0, 0.0, 2.8284271247461903]
133 - ```
134 - * diagonalization (2×2 or 3×3 matrix)
135 - ```elixir
136 - MatrixOperation.diagonalization([[1, 3], [4, 2]])
137 - [[5.0, 0], [0, -2.0]]
129 + MatrixOperation.singular_value([[1, 2, 3, 1], [2, 4, 1, 5], [3, 3, 10, 8]], 100)
130 + [14.912172620559879, 4.236463407782015, 1.6369134152873956, 0.0]
138 131 ```
139 132 * Jordan normal form (2×2 or 3×3 matrix)
140 133 ```elixir
  @@ -180,6 +173,16 @@ MatrixOperation.svd([[1, 0, 0], [0, 1, 1]], 100)
180 173 ]
181 174 ]
182 175 ```
176 + * Calculate eigenvalue by QR decomposition
177 + ```elixir
178 + MatrixOperation.eigenvalue([[6, 1, 1, 1], [1, 7, 1, 1], [1, 1, 8, 1], [1, 1, 1, 9]], 100)
179 + [10.803886359051251, 7.507748705362773, 6.39227529027387, 5.296089645312106]
180 + ```
181 + * diagonalization by QR decomposition
182 + ```elixir
183 + MatrixOperation.diagonalization([[1, 3], [4, 2]], 100)
184 + [[5.000000000000018, 0], [0, -1.999999999999997]]
185 + ```
183 186 The second argument (ex. 100) is max iterate number.
184 187 * Frobenius norm
185 188 ```elixir
  @@ -194,7 +197,7 @@ MatrixOperation.one_norm([[2, 3], [1, 4], [2, 1]])
194 197 * L two norm (L_2 norm)
195 198 ```elixir
196 199 MatrixOperation.two_norm([[2, 3], [1, 4], [2, 1]])
197 - 5.674983803488142
200 + 5.674983803488139
198 201 ```
199 202 * Max norm
200 203 ```elixir
  @@ -206,7 +209,7 @@ MatrixOperation.max_norm([[2, 3], [1, 4], [2, 1]])
206 209 MatrixOperation.rank([[2, 3, 4, 2], [1, 4, 2, 3], [2, 1, 4, 4]], 100)
207 210 3
208 211 ```
209 - The second argument (ex. 100) is max iterate number for Jacobi method.
212 + The second argument (ex. 100) is max iterate number for QR decomposition.
210 213 * Variance covariance matrix
211 214 ```elixir
212 215 MatrixOperation.variance_covariance_matrix([[40, 80], [80, 90], [90, 100]])
  @@ -11,4 +11,4 @@
11 11 <<"https://github.com/kenken-neko/elixir-matrix-operation">>}]}.
12 12 {<<"name">>,<<"matrix_operation">>}.
13 13 {<<"requirements">>,[]}.
14 - {<<"version">>,<<"0.3.5">>}.
14 + {<<"version">>,<<"0.3.6">>}.
  @@ -9,21 +9,21 @@ defmodule MatrixOperation do
9 9 iex> MatrixOperation.row_column_matrix([[3, 2, 3], [2, 1, 2]])
10 10 [2, 3]
11 11 """
12 - def row_column_matrix(a) when is_list(hd(a)) do
13 - columns_number = Enum.map(a, &row_column_matrix_sub(&1, 0))
14 - max_number = Enum.max(columns_number)
15 - if(max_number == Enum.min(columns_number), do: [length(a), max_number], else: nil)
12 + def row_column_matrix(matrix) when is_list(hd(matrix)) do
13 + col_num = Enum.map(matrix, &row_column_matrix_sub(&1, 0))
14 + max_num = Enum.max(col_num)
15 + if(max_num == Enum.min(col_num), do: [length(matrix), max_num], else: nil)
16 16 end
17 17
18 - def row_column_matrix(_) do
18 + def row_column_matrix(_matrix) do
19 19 nil
20 20 end
21 21
22 - defp row_column_matrix_sub(row_a, i) when i != length(row_a) do
23 - if(is_number(Enum.at(row_a, i)), do: row_column_matrix_sub(row_a, i + 1), else: nil)
22 + defp row_column_matrix_sub(row, i) when i != length(row) do
23 + if(is_number(Enum.at(row, i)), do: row_column_matrix_sub(row, i + 1), else: nil)
24 24 end
25 25
26 - defp row_column_matrix_sub(row_a, i) when i == length(row_a) do
26 + defp row_column_matrix_sub(row, i) when i == length(row) do
27 27 i
28 28 end
29 29
  @@ -34,8 +34,8 @@ defmodule MatrixOperation do
34 34 [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
35 35 """
36 36 def unit_matrix(n) when n > 0 and is_integer(n) do
37 - index_list = Enum.to_list(1..n)
38 - Enum.map(index_list, fn x -> Enum.map(index_list, &unit_matrix_sub(x, &1)) end)
37 + idx_list = Enum.to_list(1..n)
38 + Enum.map(idx_list, fn x -> Enum.map(idx_list, &unit_matrix_sub(x, &1)) end)
39 39 end
40 40
41 41 defp unit_matrix_sub(i, j) when i == j do
  @@ -55,11 +55,11 @@ defmodule MatrixOperation do
55 55 [[1, 1], [1, 1], [1, 1]]
56 56 """
57 57 def even_matrix(m, n, s) when m > 0 and n > 0 and is_number(s) do
58 - Enum.to_list(1..m) |>
59 - Enum.map(fn _ -> Enum.map(Enum.to_list(1..n), & &1 * 0 + s) end)
58 + Enum.to_list(1..m)
59 + |> Enum.map(fn _ -> Enum.map(Enum.to_list(1..n), & &1 * 0 + s) end)
60 60 end
61 61
62 - def even_matrix(_, _, _) do
62 + def even_matrix(_m, _n, _s) do
63 63 nil
64 64 end
65 65
  @@ -69,10 +69,10 @@ defmodule MatrixOperation do
69 69 iex> MatrixOperation.get_one_element([[1, 2, 3], [4, 5, 6], [7, 8, 9] ], [1, 1])
70 70 1
71 71 """
72 - def get_one_element(matrix, [row_index, column_index]) do
72 + def get_one_element(matrix, [row_idx, col_idx]) do
73 73 matrix
74 - |> Enum.at(row_index - 1)
75 - |> Enum.at(column_index - 1)
74 + |> Enum.at(row_idx - 1)
75 + |> Enum.at(col_idx - 1)
76 76 end
77 77
78 78 @doc """
  @@ -81,9 +81,9 @@ defmodule MatrixOperation do
81 81 iex> MatrixOperation.get_one_row([[1, 2, 3], [4, 5, 6], [7, 8, 9] ], 1)
82 82 [1, 2, 3]
83 83 """
84 - def get_one_row(matrix, row_index) do
84 + def get_one_row(matrix, row_idx) do
85 85 matrix
86 - |> Enum.at(row_index - 1)
86 + |> Enum.at(row_idx - 1)
87 87 end
88 88
89 89 @doc """
  @@ -92,10 +92,10 @@ defmodule MatrixOperation do
92 92 iex> MatrixOperation.get_one_column([[1, 2, 3], [4, 5, 6], [7, 8, 9] ], 1)
93 93 [1, 4, 7]
94 94 """
95 - def get_one_column(matrix, column_index) do
95 + def get_one_column(matrix, col_idx) do
96 96 matrix
97 - |> transpose
98 - |> Enum.at(column_index - 1)
97 + |> transpose()
98 + |> Enum.at(col_idx - 1)
99 99 end
100 100
101 101 @doc """
  @@ -104,11 +104,11 @@ defmodule MatrixOperation do
104 104 iex> MatrixOperation.delete_one_row([[1, 2, 3], [4, 5, 6], [7, 8, 9]], 3)
105 105 [[1, 2, 3], [4, 5, 6]]
106 106 """
107 - def delete_one_row(matrix, delete_index) do
107 + def delete_one_row(matrix, del_idx) do
108 108 matrix
109 109 |> Enum.with_index()
110 - |> Enum.reject(fn {_, i} -> i == delete_index - 1 end)
111 - |> Enum.map(fn {x, _} -> x end)
110 + |> Enum.reject(fn {_x, idx} -> idx == del_idx - 1 end)
111 + |> Enum.map(fn {x, _idx} -> x end)
112 112 end
113 113
114 114 @doc """
  @@ -117,13 +117,13 @@ defmodule MatrixOperation do
117 117 iex> MatrixOperation.delete_one_column([[1, 2, 3], [4, 5, 6], [7, 8, 9]], 2)
118 118 [[1, 3], [4, 6], [7, 9]]
119 119 """
120 - def delete_one_column(matrix, delete_index) do
120 + def delete_one_column(matrix, del_idx) do
121 121 matrix
122 - |> transpose
122 + |> transpose()
123 123 |> Enum.with_index()
124 - |> Enum.reject(fn {_, i} -> i == delete_index - 1 end)
125 - |> Enum.map(fn {x, _} -> x end)
126 - |> transpose
124 + |> Enum.reject(fn {_x, idx} -> idx == del_idx - 1 end)
125 + |> Enum.map(fn {x, _idx} -> x end)
126 + |> transpose()
127 127 end
128 128
129 129 @doc """
  @@ -132,10 +132,10 @@ defmodule MatrixOperation do
132 132 iex> MatrixOperation.exchange_one_row([[1, 2, 3], [4, 5, 6], [7, 8, 9]], 3, [1, 1, 1])
133 133 [[1, 2, 3], [4, 5, 6], [1, 1, 1]]
134 134 """
135 - def exchange_one_row(matrix, exchange_index, exchange_list) do
135 + def exchange_one_row(matrix, exchange_idx, exchange_list) do
136 136 matrix
137 137 |> Enum.with_index()
138 - |> Enum.map(fn {x, i} -> if(i == exchange_index - 1, do: exchange_list, else: x) end)
138 + |> Enum.map(fn {x, idx} -> if(idx == exchange_idx - 1, do: exchange_list, else: x) end)
139 139 end
140 140
141 141 @doc """
  @@ -144,27 +144,12 @@ defmodule MatrixOperation do
144 144 iex> MatrixOperation.exchange_one_column([[1, 2, 3], [4, 5, 6], [7, 8, 9]], 2, [1, 1, 1])
145 145 [[1, 1, 3], [4, 1, 6], [7, 1, 9]]
146 146 """
147 - def exchange_one_column(matrix, exchange_index, exchange_list) do
147 + def exchange_one_column(matrix, exchange_idx, exchange_list) do
148 148 matrix
149 149 |> transpose
150 150 |> Enum.with_index()
151 - |> Enum.map(fn {x, i} -> if(i == exchange_index - 1, do: exchange_list, else: x) end)
152 - |> transpose
153 - end
154 -
155 - # Count finite values
156 - defp count_finite_values(a) when is_list(a) do
157 - a
158 - |> Enum.map(&count_finite_values(&1))
159 - |> Enum.sum
160 - end
161 -
162 - defp count_finite_values(a) when is_number(a) and a == 0 do
163 - 0
164 - end
165 -
166 - defp count_finite_values(a) when is_number(a) do
167 - 1
151 + |> Enum.map(fn {x, idx} -> if(idx == exchange_idx - 1, do: exchange_list, else: x) end)
152 + |> transpose()
168 153 end
169 154
170 155 @doc """
  @@ -173,8 +158,8 @@ defmodule MatrixOperation do
173 158 iex> MatrixOperation.transpose([[1.0, 2.0], [3.0, 4.0]])
174 159 [[1.0, 3.0], [2.0, 4.0]]
175 160 """
176 - def transpose(a) do
177 - Enum.zip(a)
161 + def transpose(matrix) do
162 + Enum.zip(matrix)
178 163 |> Enum.map(&Tuple.to_list(&1))
179 164 end
180 165
  @@ -184,20 +169,20 @@ defmodule MatrixOperation do
184 169 iex> MatrixOperation.trace([[1.0, 2.0], [3.0, 4.0]])
185 170 5.0
186 171 """
187 - def trace(a) do
188 - [row, column] = row_column_matrix(a)
189 - a_index = add_index(a)
172 + def trace(matrix) do
173 + [row_num, col_num] = row_column_matrix(matrix)
174 + matrix_with_idx = add_index(matrix)
190 175
191 - Enum.map(a_index, &trace_sub(&1, row, column))
176 + Enum.map(matrix_with_idx, &trace_sub(&1, row_num, col_num))
192 177 |> Enum.sum()
193 178 end
194 179
195 - defp trace_sub(_, row, column) when row != column do
180 + defp trace_sub(_, row_num, col_num) when row_num != col_num do
196 181 nil
197 182 end
198 183
199 - defp trace_sub([index, row_list], _, _) do
200 - Enum.at(row_list, index - 1)
184 + defp trace_sub([idx, row], _row_num, _col_num) do
185 + Enum.at(row, idx - 1)
201 186 end
202 187
203 188 @doc """
  @@ -210,20 +195,13 @@ defmodule MatrixOperation do
210 195 iex> MatrixOperation.determinant([ [3,1,1,2,1], [5,1,3,4,1], [2,0,1,0,1], [1,3,2,1,1], [1,1,1,1,1] ])
211 196 -14
212 197 """
213 - def determinant(a) do
214 - determinant_sub(1, a)
215 - end
216 -
217 - # minor_matrix
218 - defp minor_matrix(a_with_index, row) do
219 - (a_with_index -- [row])
220 - |> Enum.map(&Enum.at(&1, 1))
221 - |> Enum.map(&Enum.drop(&1, 1))
198 + def determinant(matrix) do
199 + determinant_sub(1, matrix)
222 200 end
223 201
224 202 # 1×1 matrix
225 - defp determinant_sub(_, a) when length(a) == 1 do
226 - Enum.at(a, 0)
203 + defp determinant_sub(_, matrix) when length(matrix) == 1 do
204 + Enum.at(matrix, 0)
227 205 |> Enum.at(0)
228 206 end
229 207
  @@ -233,23 +211,29 @@ defmodule MatrixOperation do
233 211 end
234 212
235 213 # 3×3 or over matrix
236 - defp determinant_sub(co, a) do
237 - a_with_index = add_index(a)
214 + defp determinant_sub(co, matrix) do
215 + matrix_with_idx = add_index(matrix)
238 216
239 217 Enum.map(
240 - a_with_index,
218 + matrix_with_idx,
241 219 &determinant_sub(
242 220 (-1 + 2 * rem(hd(&1), 2)) * co * hd(Enum.at(&1, 1)),
243 - minor_matrix(a_with_index, &1)
221 + minor_matrix(matrix_with_idx, &1)
244 222 )
245 223 )
246 224 |> Enum.sum()
247 225 end
248 226
227 + defp minor_matrix(matrix_with_idx, row) do
228 + (matrix_with_idx -- [row])
229 + |> Enum.map(&Enum.at(&1, 1))
230 + |> Enum.map(&Enum.drop(&1, 1))
231 + end
232 +
249 233 # add index
250 - defp add_index(a) do
234 + defp add_index(matrix) do
251 235 Stream.iterate(1, &(&1 + 1))
252 - |> Enum.zip(a)
236 + |> Enum.zip(matrix)
253 237 |> Enum.map(&(&1 |> Tuple.to_list()))
254 238 end
255 239
  @@ -261,36 +245,36 @@ defmodule MatrixOperation do
261 245 iex> MatrixOperation.cramer([[0, -2, 1], [-1, 1, -4], [3, 3, 1]], [[3], [-7], [4]], 1)
262 246 2.0
263 247 """
264 - def cramer(a, vertical_vec, select_index) do
248 + def cramer(matrix, vertical_vec, select_idx) do
265 249 [t] = transpose(vertical_vec)
266 - det_a = determinant(a)
267 - cramer_sub(a, t, select_index - 1, det_a)
250 + det = determinant(matrix)
251 + cramer_sub(matrix, t, select_idx - 1, det)
268 252 end
269 253
270 254 defp cramer_sub(_, _, _, nil), do: nil
271 255 defp cramer_sub(_, _, _, 0), do: nil
272 256
273 - defp cramer_sub(a, t, select_index, det_a) do
274 - rep_det_a = transpose(a) |> replace_element_in_list(select_index, t, 0, []) |> determinant
275 - rep_det_a / det_a
257 + defp cramer_sub(a, t, select_idx, det) do
258 + rep_det = transpose(a) |> replace_element_in_list(select_idx, t, 0, []) |> determinant
259 + rep_det / det
276 260 end
277 261
278 262 defp replace_element_in_list(list, i, replace_element, i, output) when i < length(list) do
279 263 replace_element_in_list(list, i, replace_element, i + 1, output ++ [replace_element])
280 264 end
281 265
282 - defp replace_element_in_list(list, select_index, replace_element, i, output)
266 + defp replace_element_in_list(list, select_idx, replace_element, i, output)
283 267 when i < length(list) do
284 268 replace_element_in_list(
285 269 list,
286 - select_index,
270 + select_idx,
287 271 replace_element,
288 272 i + 1,
289 273 output ++ [Enum.at(list, i)]
290 274 )
291 275 end
292 276
293 - defp replace_element_in_list(list, _select_index, _replace_element, i, output)
277 + defp replace_element_in_list(list, _select_idx, _replace_element, i, output)
294 278 when i == length(list),
295 279 do: output
296 280
  @@ -302,22 +286,22 @@ defmodule MatrixOperation do
302 286 iex> MatrixOperation.linear_equations_cramer([[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[1], [0], [0]])
303 287 [1.0, 0.0, 0.0]
304 288 """
305 - def linear_equations_cramer(a, vertical_vec) do
289 + def linear_equations_cramer(matrix, vertical_vec) do
306 290 # check the setupufficient condition
307 - if determinant(a) == 0 do
291 + if determinant(matrix) == 0 do
308 292 nil
309 293 else
310 294 [t] = transpose(vertical_vec)
311 - linear_equations_cramer_sub(a, t, 0, [])
295 + linear_equations_cramer_sub(matrix, t, 0, [])
312 296 end
313 297 end
314 298
315 - defp linear_equations_cramer_sub(a, t, i, output) when i < length(a) do
299 + defp linear_equations_cramer_sub(matrix, t, i, output) when i < length(matrix) do
316 300 vertical_vec = transpose([t])
317 - linear_equations_cramer_sub(a, t, i + 1, output ++ [cramer(a, vertical_vec, i + 1)])
301 + linear_equations_cramer_sub(matrix, t, i + 1, output ++ [cramer(matrix, vertical_vec, i + 1)])
318 302 end
319 303
320 - defp linear_equations_cramer_sub(a, _t, i, output) when i == length(a) do
304 + defp linear_equations_cramer_sub(matrix, _t, i, output) when i == length(matrix) do
321 305 output
322 306 end
323 307
  @@ -327,8 +311,9 @@ defmodule MatrixOperation do
327 311 iex> MatrixOperation.leading_principal_minor([[1, 3, 2], [2, 5, 1], [3, 4, 5]], 2)
328 312 [[1, 3], [2, 5]]
329 313 """
330 - def leading_principal_minor(a, k) do
331 - Enum.slice(a, 0, k)
314 + def leading_principal_minor(matrix, k) do
315 + matrix
316 + |> Enum.slice(0, k)
332 317 |> Enum.map(& Enum.slice(&1, 0, k))
333 318 end
334 319
  @@ -341,56 +326,57 @@ defmodule MatrixOperation do
341 326 [[1, 1, 0, 3], [0, -1.0, -1.0, -5.0], [0, 0, 3.0, 13.0], [0, 0, 0, -13.0]]
342 327 ]
343 328 """
344 - def lu_decomposition(a) do
345 - row_column = row_column_matrix(a)
329 + def lu_decomposition(matrix) do
330 + [row_num, col_num] = row_column_matrix(matrix)
346 331 # check the setupufficient condition
347 - check_number = lu_decomposition_check(a, row_column)
348 - if(check_number == 0, do: nil, else: lu_decomposition_sub(a, 0, length(a), [], []))
332 + check_num = lu_decomposition_check(matrix, row_num, col_num)
333 + if(check_num == 0, do: nil, else: lu_decomposition_sub(matrix, 0, length(matrix), [], []))
349 334 end
350 335
351 - defp lu_decomposition_check(_, [row_num, column_num]) when row_num != column_num do
336 + defp lu_decomposition_check(_matrix, row_num, col_num) when row_num != col_num do
352 337 nil
353 338 end
354 339
355 - defp lu_decomposition_check(a, [row_num, _]) do
340 + defp lu_decomposition_check(matrix, row_num, _col_num) do
356 341 Enum.to_list(1..row_num)
357 - |> Enum.map(& leading_principal_minor(a, &1) |> determinant)
342 + |> Enum.map(& leading_principal_minor(matrix, &1) |> determinant)
358 343 |> Enum.reduce(fn x, acc -> x * acc end)
359 344 end
360 345
361 - defp lu_decomposition_sub(a, k, len_a, _, _) when k == 0 do
362 - u_matrix = even_matrix(len_a, len_a, 0)
363 - |> exchange_one_row(1, hd(a))
346 + defp lu_decomposition_sub(matrix, k, matrix_len, _l_matrix, _u_matrix) when k == 0 do
347 + u_matrix = even_matrix(matrix_len, matrix_len, 0)
348 + |> exchange_one_row(1, hd(matrix))
364 349 inverce_u11 = 1.0 / hd(hd(u_matrix))
365 - a_factor = transpose(a)
366 - |> get_one_row(1)
367 - |> Enum.slice(1, len_a)
368 - l_row = [1] ++ hd(const_multiple(inverce_u11, [a_factor]))
369 - l_matrix = even_matrix(len_a, len_a, 0)
350 + factor = matrix
351 + |> transpose()
352 + |> get_one_row(1)
353 + |> Enum.slice(1, matrix_len)
354 + l_row = [1] ++ hd(const_multiple(inverce_u11, [factor]))
355 + l_matrix = even_matrix(matrix_len, matrix_len, 0)
370 356 |> exchange_one_row(1, l_row)
371 - lu_decomposition_sub(a, k + 1, len_a, l_matrix, u_matrix)
357 + lu_decomposition_sub(matrix, k + 1, matrix_len, l_matrix, u_matrix)
372 358 end
373 359
374 - defp lu_decomposition_sub(a, k, len_a, l_matrix, u_matrix) when k != len_a do
375 - a_t = transpose(a)
376 - u_solve = u_cal(a, k, len_a, l_matrix, u_matrix)
360 + defp lu_decomposition_sub(matrix, k, matrix_len, l_matrix, u_matrix) when k != matrix_len do
361 + t_matrix = transpose(matrix)
362 + u_solve = u_cal(matrix, k, matrix_len, l_matrix, u_matrix)
377 363 u_matrix_2 = exchange_one_row(u_matrix, k + 1, u_solve)
378 - l_solve = l_cal(a_t, k, len_a, l_matrix, u_matrix_2)
364 + l_solve = l_cal(t_matrix, k, matrix_len, l_matrix, u_matrix_2)
379 365 l_matrix_2 = exchange_one_row(l_matrix, k + 1, l_solve)
380 - lu_decomposition_sub(a, k + 1, len_a, l_matrix_2, u_matrix_2)
366 + lu_decomposition_sub(matrix, k + 1, matrix_len, l_matrix_2, u_matrix_2)
381 367 end
382 368
383 - defp lu_decomposition_sub(_, _, _, l_matrix, u_matrix) do
369 + defp lu_decomposition_sub(_matrix, _k, _matrix_len, l_matrix, u_matrix) do
384 370 [transpose(l_matrix), u_matrix]
385 371 end
386 372
387 - defp l_cal(a_t, k, len_a, l_matrix, u_matrix) do
388 - a_factor = Enum.at(a_t, k) |> Enum.slice(k + 1, len_a)
373 + defp l_cal(t_matrix, k, matrix_len, l_matrix, u_matrix) do
374 + factor = Enum.at(t_matrix, k) |> Enum.slice(k + 1, matrix_len)
389 375 u_extract = transpose(u_matrix) |> Enum.at(k)
390 376 l_row = transpose(l_matrix)
391 - |> Enum.slice(k + 1, len_a)
377 + |> Enum.slice(k + 1, matrix_len)
392 378 |> Enum.map(& inner_product(&1, u_extract))
393 - |> Enum.zip(a_factor)
379 + |> Enum.zip(factor)
394 380 |> Enum.map(fn {x, y} -> y - x end)
395 381
396 382 inverce_uii = 1.0 / Enum.at(Enum.at(u_matrix, k), k)
  @@ -399,13 +385,13 @@ defmodule MatrixOperation do
399 385 |> add_zero_element(0, k)
400 386 end
401 387
402 - defp u_cal(a, k, len_a, l_matrix, u_matrix) do
403 - a_factor = Enum.at(a, k) |> Enum.slice(k, len_a)
388 + defp u_cal(matrix, k, matrix_len, l_matrix, u_matrix) do
389 + factor = Enum.at(matrix, k) |> Enum.slice(k, matrix_len)
404 390 l_extract = transpose(l_matrix) |> Enum.at(k)
405 391 transpose(u_matrix)
406 - |> Enum.slice(k, len_a)
392 + |> Enum.slice(k, matrix_len)
407 393 |> Enum.map(& inner_product(&1, l_extract))
408 - |> Enum.zip(a_factor)
394 + |> Enum.zip(factor)
409 395 |> Enum.map(fn {x, y} -> y - x end)
410 396 |> add_zero_element(0, k)
411 397 end
  @@ -414,7 +400,7 @@ defmodule MatrixOperation do
414 400 add_zero_element([0] ++ list, init + 1, fin)
415 401 end
416 402
417 - defp add_zero_element(list, _, _) do
403 + defp add_zero_element(list, _init, _fin) do
418 404 list
419 405 end
420 406
  @@ -426,24 +412,24 @@ defmodule MatrixOperation do
426 412 iex> MatrixOperation.linear_equations_direct([[4, 1, 1], [1, 3, 1], [2, 1, 5]], [[9], [10], [19]])
427 413 [1.0, 2.0, 3.0]
428 414 """
429 - def linear_equations_direct(a, vertical_vec) do
415 + def linear_equations_direct(matrix, vertical_vec) do
430 416 # check the setupufficient condition
431 - if determinant(a) == 0 do
417 + if determinant(matrix) == 0 do
432 418 nil
433 419 else
434 420 [t] = transpose(vertical_vec)
435 - lu_decomposition_sub(a, t)
421 + linear_equations_direct_sub(matrix, t)
436 422 end
437 423 end
438 424
439 - defp lu_decomposition_sub(a, t) do
440 - [l_matrix, u_matrix] = lu_decomposition(a)
425 + defp linear_equations_direct_sub(matrix, t) do
426 + [l_matrix, u_matrix] = lu_decomposition(matrix)
441 427 dim = length(l_matrix)
442 428 y = forward_substitution(l_matrix, t, [], 0, dim)
443 429 backward_substitution(u_matrix, y, [], dim, dim)
444 430 end
445 431
446 - defp forward_substitution(l_matrix, t, _, k, dim) when k == 0 do
432 + defp forward_substitution(l_matrix, t, _y, k, dim) when k == 0 do
447 433 forward_substitution(l_matrix, t, [hd(t)], k + 1, dim)
448 434 end
449 435
  @@ -455,11 +441,11 @@ defmodule MatrixOperation do
455 441 forward_substitution(l_matrix, t, y ++ [t_ly], k + 1, dim)
456 442 end
457 443
458 - defp forward_substitution(_, _, y, k, dim) when k == dim do
444 + defp forward_substitution(_l_matrix, _t, y, k, dim) when k == dim do
459 445 y
460 446 end
461 447
462 - defp backward_substitution(u_matrix, y, _, k, dim) when k == dim do
448 + defp backward_substitution(u_matrix, y, _b, k, dim) when k == dim do
463 449 dim_1 = dim - 1
464 450 y_n = Enum.at(y, dim_1)
465 451 u_nn = Enum.at(Enum.at(u_matrix, dim_1), dim_1)
  @@ -487,12 +473,12 @@ defmodule MatrixOperation do
487 473 iex> MatrixOperation.const_multiple(2, [[1, 2, 3], [2, 2, 2], [3, 8, 9]])
488 474 [[2, 4, 6], [4, 4, 4], [6, 16, 18]]
489 475 """
490 - def const_multiple(const, a) when is_number(a) do
491 - const * a
476 + def const_multiple(const, x) when is_number(x) do
477 + const * x
492 478 end
493 479
494 - def const_multiple(const, a) when is_list(a) do
495 - Enum.map(a, &const_multiple(const, &1))
480 + def const_multiple(const, x) when is_list(x) do
481 + Enum.map(x, &const_multiple(const, &1))
496 482 end
497 483
498 484 @doc """
  @@ -501,12 +487,12 @@ defmodule MatrixOperation do
501 487 iex> MatrixOperation.const_addition(1, [1.0, 2.0, 3.0])
502 488 [2.0, 3.0, 4.0]
503 489 """
504 - def const_addition(const, a) when is_number(a) do
505 - const + a
490 + def const_addition(const, x) when is_number(x) do
491 + const + x
506 492 end
507 493
508 - def const_addition(const, a) when is_list(a) do
509 - Enum.map(a, &const_addition(const, &1))
494 + def const_addition(const, x) when is_list(x) do
495 + Enum.map(x, &const_addition(const, &1))
510 496 end
511 497
512 498 @doc """
  @@ -515,42 +501,41 @@ defmodule MatrixOperation do
515 501 iex> MatrixOperation.inverse_matrix([[1, 1, -1], [-2, -1, 1], [-1, -2, 1]])
516 502 [[-1.0, -1.0, 0.0], [-1.0, 0.0, -1.0], [-3.0, -1.0, -1.0]]
517 503 """
518 - def inverse_matrix(a) when is_list(hd(a)) do
519 - det_a = determinant(a)
504 + def inverse_matrix(matrix) when is_list(hd(matrix)) do
505 + det = determinant(matrix)
520 506
521 - create_index_matrix(a)
522 - |> Enum.map(&map_index_row(a, det_a, &1))
523 - |> transpose
507 + create_index_matrix(matrix)
508 + |> Enum.map(&map_index_row(matrix, det, &1))
509 + |> transpose()
524 510 end
525 511
526 512 def inverse_matrix(_) do
527 513 nil
528 514 end
529 515
530 - defp create_index_matrix(a) do
531 - index_list = Enum.to_list(1..length(a))
532 - Enum.map(index_list, fn x -> Enum.map(index_list, &[x, &1]) end)
516 + defp create_index_matrix(matrix) do
517 + idx_list = Enum.to_list(1..length(matrix))
518 + Enum.map(idx_list, fn x -> Enum.map(idx_list, &[x, &1]) end)
533 519 end
534 520
535 - defp map_index_row(_, 0, _) do
521 + defp map_index_row(_matrix, det, _row) when det == 0 do
536 522 nil
537 523 end
538 524
539 - defp map_index_row(a, det_a, row) do
540 - Enum.map(row, &minor_matrix(a, det_a, &1))
525 + defp map_index_row(matrix, det, row) do
526 + Enum.map(row, &minor_matrix(matrix, det, &1))
541 527 end
542 528
543 - # minor_matrix
544 - defp minor_matrix(a, det_a, [row_number, column_number]) do
545 - det_temp_a =
546 - delete_one_row(a, row_number)
529 + defp minor_matrix(matrix, det, [row_num, col_num]) do
530 + det_temp_matrix =
531 + delete_one_row(matrix, row_num)
547 532 |> transpose
548 - |> delete_one_row(column_number)
533 + |> delete_one_row(col_num)
549 534 |> determinant
550 535
551 - if(rem(row_number + column_number, 2) == 0,
552 - do: det_temp_a / det_a,
553 - else: -1 * det_temp_a / det_a
536 + if(rem(row_num + col_num, 2) == 0,
537 + do: det_temp_matrix / det,
538 + else: -1 * det_temp_matrix / det
554 539 )
555 540 end
556 541
  @@ -565,9 +550,9 @@ defmodule MatrixOperation do
565 550 end
566 551
567 552 defp check_product(a, b) do
568 - column_number_a = row_column_matrix(a) |> Enum.at(1)
569 - row_number_b = row_column_matrix(b) |> Enum.at(0)
570 - if(column_number_a == row_number_b, do: product_sub(a, b), else: nil)
553 + col_num_a = row_column_matrix(a) |> Enum.at(1)
554 + row_num_b = row_column_matrix(b) |> Enum.at(0)
555 + if(col_num_a == row_num_b, do: product_sub(a, b), else: nil)
571 556 end
572 557
573 558 defp product_sub(a, b) do
  @@ -577,8 +562,8 @@ defmodule MatrixOperation do
577 562 end)
578 563 end
579 564
580 - defp inner_product(row_a, column_b) do
581 - Enum.zip(row_a, column_b)
565 + defp inner_product(row_a, col_b) do
566 + Enum.zip(row_a, col_b)
582 567 |> Enum.map(&Tuple.to_list(&1))
583 568 |> Enum.map(&Enum.reduce(&1, fn x, acc -> x * acc end))
584 569 |> Enum.sum()
  @@ -595,9 +580,9 @@ defmodule MatrixOperation do
595 580 end
596 581
597 582 defp check_add(a, b) do
598 - row_column_a = row_column_matrix(a)
599 - row_column_b = row_column_matrix(b)
600 - if(row_column_a == row_column_b, do: add_sub(a, b), else: nil)
583 + row_col_a = row_column_matrix(a)
584 + row_col_b = row_column_matrix(b)
585 + if(row_col_a == row_col_b, do: add_sub(a, b), else: nil)
601 586 end
602 587
603 588 defp add_sub(a, b) do
  @@ -620,9 +605,9 @@ defmodule MatrixOperation do
620 605 end
621 606
622 607 defp check_subtract(a, b) do
623 - row_column_a = row_column_matrix(a)
624 - row_column_b = row_column_matrix(b)
625 - if(row_column_a == row_column_b, do: subtract_sub(a, b), else: nil)
608 + row_col_a = row_column_matrix(a)
609 + row_col_b = row_column_matrix(b)
610 + if(row_col_a == row_col_b, do: subtract_sub(a, b), else: nil)
626 611 end
627 612
628 613 defp subtract_sub(a, b) do
  @@ -674,8 +659,8 @@ defmodule MatrixOperation do
674 659 iex> MatrixOperation.hadamard_power([[3, 2, 3], [2, 1, 2]], 2)
675 660 [[9.0, 4.0, 9.0], [4.0, 1.0, 4.0]]
676 661 """
677 - def hadamard_power(a, n) do
678 - Enum.map(a, &Enum.map(&1, fn x -> :math.pow(x, n) end))
662 + def hadamard_power(matrix, n) do
663 + Enum.map(matrix, &Enum.map(&1, fn x -> :math.pow(x, n) end))
679 664 end
680 665
681 666 @doc """
  @@ -704,24 +689,24 @@ defmodule MatrixOperation do
704 689 end
705 690
706 691 @doc """
707 - eigenvalue by the direct method [R^2×R^2/R^3×R^3 matrix]
692 + eigenvalue by algebra method [R^2×R^2/R^3×R^3 matrix]
708 693 ## Examples
709 - iex> MatrixOperation.eigenvalue_direct([[3, 1], [2, 2]])
694 + iex> MatrixOperation.eigenvalue_algebra([[3, 1], [2, 2]])
710 695 [4.0, 1.0]
711 - iex> MatrixOperation.eigenvalue_direct([[6, -3], [4, -1]])
696 + iex> MatrixOperation.eigenvalue_algebra([[6, -3], [4, -1]])
712 697 [3.0, 2.0]
713 - iex> MatrixOperation.eigenvalue_direct([[1, 1, 1], [1, 2, 1], [1, 2, 3]])
698 + iex> MatrixOperation.eigenvalue_algebra([[1, 1, 1], [1, 2, 1], [1, 2, 3]])
714 699 [4.561552806429505, 0.43844714673139706, 1.0000000468390973]
715 - iex> MatrixOperation.eigenvalue_direct([[2, 1, -1], [1, 1, 0], [-1, 0, 1]])
700 + iex> MatrixOperation.eigenvalue_algebra([[2, 1, -1], [1, 1, 0], [-1, 0, 1]])
716 701 [3.0000000027003626, 0, 0.9999999918989121]
717 702 """
718 703 # 2×2 algebra method
719 - def eigenvalue_direct([[a11, a12], [a21, a22]]) do
704 + def eigenvalue_algebra([[a11, a12], [a21, a22]]) do
720 705 quadratic_formula(1, -a11 - a22, a11 * a22 - a12 * a21)
721 706 end
722 707
723 708 # 3×3 algebratic method
724 - def eigenvalue_direct([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]) do
709 + def eigenvalue_algebra([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]) do
725 710 a = -1
726 711 b = a11 + a22 + a33
727 712 c = a21 * a12 + a13 * a31 + a32 * a23 - a11 * a22 - a11 * a33 - a22 * a33
  @@ -733,7 +718,7 @@ defmodule MatrixOperation do
733 718 if(dis > 0, do: cubic_formula(a, b, c, d), else: nil)
734 719 end
735 720
736 - def eigenvalue_direct(_a) do
721 + def eigenvalue_algebra(_a) do
737 722 "2×2 or 3×3 matrix only"
738 723 end
739 724
  @@ -851,38 +836,38 @@ defmodule MatrixOperation do
851 836 end
852 837
853 838 @doc """
854 - Matrix diagonalization by the direct method [R^2×R^2/R^3×R^3 matrix]
839 + Matrix diagonalization by algebra method [R^2×R^2/R^3×R^3 matrix]
855 840 #### Examples
856 - iex> MatrixOperation.diagonalization_direct([[1, 3], [4, 2]])
841 + iex> MatrixOperation.diagonalization_algebra([[1, 3], [4, 2]])
857 842 [[5.0, 0], [0, -2.0]]
858 - iex> MatrixOperation.diagonalization_direct([[2, 1, -1], [1, 1, 0], [-1, 0, 1]])
843 + iex> MatrixOperation.diagonalization_algebra([[2, 1, -1], [1, 1, 0], [-1, 0, 1]])
859 844 [[3.0000000027003626, 0, 0], [0, 0, 0], [0, 0, 0.9999999918989121]]
860 845 """
861 - def diagonalization_direct(a) do
862 - eigenvalue_direct(a)
863 - |> diagonalization_direct_condition
846 + def diagonalization_algebra(matrix) do
847 + eigenvalue_algebra(matrix)
848 + |> diagonalization_algebra_condition()
864 849 end
865 850
866 - defp diagonalization_direct_condition(a) when a == nil do
851 + defp diagonalization_algebra_condition(matrix) when matrix == nil do
867 852 nil
868 853 end
869 854
870 - defp diagonalization_direct_condition(a) do
871 - a
872 - |> Enum.with_index
873 - |> Enum.map(& diagonalization_direct_sub(&1, length(a), 0, []))
855 + defp diagonalization_algebra_condition(matrix) do
856 + matrix
857 + |> Enum.with_index()
858 + |> Enum.map(& diagonalization_algebra_sub(&1, length(matrix), 0, []))
874 859 end
875 860
876 - defp diagonalization_direct_sub(_, dim, i, row) when i + 1 > dim do
861 + defp diagonalization_algebra_sub(_, dim, i, row) when i + 1 > dim do
877 862 row
878 863 end
879 864
880 - defp diagonalization_direct_sub({ev, index}, dim, i, row) when i != index do
881 - diagonalization_direct_sub({ev, index}, dim, i + 1, row ++ [0])
865 + defp diagonalization_algebra_sub({ev, index}, dim, i, row) when i != index do
866 + diagonalization_algebra_sub({ev, index}, dim, i + 1, row ++ [0])
882 867 end
883 868
884 - defp diagonalization_direct_sub({ev, index}, dim, i, row) when i == index do
885 - diagonalization_direct_sub({ev, index}, dim, i + 1, row ++ [ev])
869 + defp diagonalization_algebra_sub({ev, index}, dim, i, row) when i == index do
870 + diagonalization_algebra_sub({ev, index}, dim, i + 1, row ++ [ev])
886 871 end
887 872
888 873 @doc """
  @@ -924,7 +909,7 @@ defmodule MatrixOperation do
924 909 end
925 910
926 911 defp jordan_R2R2(b, c, m) when (b * b > 4 * c) do
927 - diagonalization_direct(m)
912 + diagonalization_algebra(m)
928 913 end
929 914
930 915 defp jordan_R2R2(b, c, m) when b * b == 4 * c do
  @@ -965,7 +950,7 @@ defmodule MatrixOperation do
965 950 defp jordan_R3R3(b, c, d, m)
966 951 when 4 * c * c * c - 27 * d * d + b * b * c * c - 18 * b * c * d -
967 952 4 * b * b * b * d > 0 do
968 - diagonalization_direct(m)
953 + diagonalization_algebra(m)
969 954 end
970 955 # Triple root
971 956 defp jordan_R3R3(b, c, d, m)
  @@ -1047,11 +1032,11 @@ defmodule MatrixOperation do
1047 1032 [1.0, -2.0, 2.0]
1048 1033 ]
1049 1034 """
1050 - def power_iteration(a, max_k) do
1051 - init_vec = random_column(length(a))
1052 - xk_pre = power_iteration_sub(a, init_vec, max_k)
1035 + def power_iteration(matrix, iter_num) do
1036 + init_vec = random_column(length(matrix))
1037 + xk_pre = power_iteration_sub(matrix, init_vec, iter_num)
1053 1038 # eigen vector
1054 - [xk_vec] = product(a, xk_pre) |> transpose
1039 + [xk_vec] = product(matrix, xk_pre) |> transpose
1055 1040 [xk_pre_vec] = transpose(xk_pre)
1056 1041 # eigen value
1057 1042 eigen_value = inner_product(xk_vec, xk_vec) / inner_product(xk_vec, xk_pre_vec)
  @@ -1067,10 +1052,10 @@ defmodule MatrixOperation do
1067 1052 nil
1068 1053 end
1069 1054
1070 - defp power_iteration_sub(a, v, max_k) do
1055 + defp power_iteration_sub(matrix, v, iter_num) do
1071 1056 # Normarization is for overflow suppression
1072 - Enum.reduce(1..max_k, v, fn _, acc ->
1073 - vp = product(a, acc)
1057 + Enum.reduce(1..iter_num, v, fn _, acc ->
1058 + vp = product(matrix, acc)
1074 1059 [vpt] = transpose(vp)
1075 1060 const_multiple(1 / :math.sqrt(inner_product(vpt, vpt)), vp)
1076 1061 end)
  @@ -1089,43 +1074,43 @@ defmodule MatrixOperation do
1089 1074 ]
1090 1075 ]
1091 1076 """
1092 - def jacobi(a, loop_num) do
1093 - [pap, p] = jacobi_loop(a, loop_num, 0, unit_matrix(length(a)))
1077 + def jacobi(matrix, iter_num) do
1078 + [pap, p] = jacobi_iteration(matrix, iter_num, 0, unit_matrix(length(matrix)))
1094 1079 p_rnd = Enum.map(p, & Enum.map(&1, fn x -> zero_approximation(x) end))
1095 1080
1096 1081 eigenvalue_list = pap
1097 - |> Enum.with_index
1082 + |> Enum.with_index()
1098 1083 |> Enum.map(& jacobi_sub4(&1))
1099 1084 |> Enum.map(& zero_approximation(&1))
1100 1085 [eigenvalue_list, p_rnd]
1101 1086 end
1102 1087
1103 - defp jacobi_loop(a, loop_num, l, p_pre) when l != loop_num do
1104 - [row_num, column_num] = row_column_matrix(a)
1105 - odts = off_diagonal_terms(a, row_num, column_num, 0, 0, [])
1088 + defp jacobi_iteration(matrix, iter_num, l, p_pre) when l != iter_num do
1089 + [row_num, col_num] = row_column_matrix(matrix)
1090 + odts = off_diagonal_terms(matrix, row_num, col_num, 0, 0, [])
1106 1091 |> Enum.map(& abs(&1))
1107 1092
1108 1093 max_odt = Enum.max(odts)
1109 1094 [max_i, max_j] = Enum.with_index(odts)
1110 1095 |> jocobi_sub(max_odt, 0)
1111 - |> jocobi_sub2(column_num, 0)
1096 + |> jocobi_sub2(col_num, 0)
1112 1097
1113 - a_ij = get_one_element(a, [max_i + 1, max_j + 1])
1114 - a_ii = get_one_element(a, [max_i + 1, max_i + 1])
1115 - a_jj = get_one_element(a, [max_j + 1, max_j + 1])
1098 + a_ij = get_one_element(matrix, [max_i + 1, max_j + 1])
1099 + a_ii = get_one_element(matrix, [max_i + 1, max_i + 1])
1100 + a_jj = get_one_element(matrix, [max_j + 1, max_j + 1])
1116 1101 phi = phi_if(a_ii - a_jj, a_ij)
1117 1102
1118 - p = jacobi_sub3(phi, column_num, max_i, max_j, 0, 0, [], [])
1103 + p = jacobi_sub3(phi, col_num, max_i, max_j, 0, 0, [], [])
1119 1104 p_pi = product(p_pre, p)
1120 1105 p
1121 - |> transpose
1122 - |> product(a)
1106 + |> transpose()
1107 + |> product(matrix)
1123 1108 |> product(p)
1124 - |> jacobi_loop(loop_num, l + 1, p_pi)
1109 + |> jacobi_iteration(iter_num, l + 1, p_pi)
1125 1110 end
1126 1111
1127 - defp jacobi_loop(a, _, _, p) do
1128 - [a, p]
1112 + defp jacobi_iteration(matrix, _, _, p) do
1113 + [matrix, p]
1129 1114 end
1130 1115
1131 1116 defp phi_if(denominator, a_ij) when denominator < 0.0000001 and a_ij > 0 do
  @@ -1140,20 +1125,20 @@ defmodule MatrixOperation do
1140 1125 atan(-2 * a_ij / denominator) * 0.5
1141 1126 end
1142 1127
1143 - defp off_diagonal_terms(m, row_num, column_num, i, j, output) when i < j and row_num >= i and column_num > j do
1144 - off_diagonal_terms(m, row_num, column_num, i, j + 1, output ++ [get_one_element(m, [i + 1, j + 1])])
1128 + defp off_diagonal_terms(m, row_num, col_num, i, j, output) when i < j and row_num >= i and col_num > j do
1129 + off_diagonal_terms(m, row_num, col_num, i, j + 1, output ++ [get_one_element(m, [i + 1, j + 1])])
1145 1130 end
1146 1131
1147 - defp off_diagonal_terms(m, row_num, column_num, i, j, output) when i < j and row_num > i and column_num == j do
1148 - off_diagonal_terms(m, row_num, column_num, i + 1, 0, output)
1132 + defp off_diagonal_terms(m, row_num, col_num, i, j, output) when i < j and row_num > i and col_num == j do
1133 + off_diagonal_terms(m, row_num, col_num, i + 1, 0, output)
1149 1134 end
1150 1135
1151 - defp off_diagonal_terms(_, row_num, column_num, i, j, output) when row_num == i and column_num == j do
1136 + defp off_diagonal_terms(_, row_num, col_num, i, j, output) when row_num == i and col_num == j do
1152 1137 output
1153 1138 end
1154 1139
1155 - defp off_diagonal_terms(m, row_num, column_num, i, j, output) do
1156 - off_diagonal_terms(m, row_num, column_num, i, j + 1, output)
1140 + defp off_diagonal_terms(m, row_num, col_num, i, j, output) do
1141 + off_diagonal_terms(m, row_num, col_num, i, j + 1, output)
1157 1142 end
1158 1143
1159 1144 defp jocobi_sub(element_idx_list, target_element, i) when hd(element_idx_list) == {target_element, i} do
  @@ -1165,45 +1150,45 @@ defmodule MatrixOperation do
1165 1150 jocobi_sub(tail, target_element, i + 1)
1166 1151 end
1167 1152
1168 - defp jocobi_sub2(idx, column_num, i) when idx < (i + 1) * column_num - ((i + 1) * (2 + i) * 0.5) do
1169 - [max_i, max_j] = [i, idx - i * (2 * column_num - i - 1) * 0.5 + i + 1]
1153 + defp jocobi_sub2(idx, col_num, i) when idx < (i + 1) * col_num - ((i + 1) * (2 + i) * 0.5) do
1154 + [max_i, max_j] = [i, idx - i * (2 * col_num - i - 1) * 0.5 + i + 1]
1170 1155 [max_i, round(max_j)]
1171 1156 end
1172 1157
1173 - defp jocobi_sub2(idx, column_num, i) do
1174 - jocobi_sub2(idx, column_num, i + 1)
1158 + defp jocobi_sub2(idx, col_num, i) do
1159 + jocobi_sub2(idx, col_num, i + 1)
1175 1160 end
1176 1161
1177 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when i == j and ( i == target_i or j == target_j) do
1178 - jacobi_sub3(phi, column_num, target_i, target_j, i, j + 1, o_row ++ [:math.cos(phi)], output)
1162 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when i == j and ( i == target_i or j == target_j) do
1163 + jacobi_sub3(phi, col_num, target_i, target_j, i, j + 1, o_row ++ [:math.cos(phi)], output)
1179 1164 end
1180 1165
1181 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when i == target_i and j == target_j and j != column_num do
1182 - jacobi_sub3(phi, column_num, target_i, target_j, i, j + 1, o_row ++ [:math.sin(phi)], output)
1166 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when i == target_i and j == target_j and j != col_num do
1167 + jacobi_sub3(phi, col_num, target_i, target_j, i, j + 1, o_row ++ [:math.sin(phi)], output)
1183 1168 end
1184 1169
1185 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when i == target_i and j == target_j and j == column_num do
1186 - jacobi_sub3(phi, column_num, target_i, target_j, i + 1, 0, [] , output ++ [o_row ++ [:math.sin(phi)]])
1170 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when i == target_i and j == target_j and j == col_num do
1171 + jacobi_sub3(phi, col_num, target_i, target_j, i + 1, 0, [] , output ++ [o_row ++ [:math.sin(phi)]])
1187 1172 end
1188 1173
1189 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when i == target_j and j == target_i do
1190 - jacobi_sub3(phi, column_num, target_i, target_j, i, j + 1, o_row ++ [:math.sin(-phi)], output)
1174 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when i == target_j and j == target_i do
1175 + jacobi_sub3(phi, col_num, target_i, target_j, i, j + 1, o_row ++ [:math.sin(-phi)], output)
1191 1176 end
1192 1177
1193 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when (i != target_i or j != target_j) and i == j and j != column_num do
1194 - jacobi_sub3(phi, column_num, target_i, target_j, i, j + 1, o_row ++ [1], output)
1178 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when (i != target_i or j != target_j) and i == j and j != col_num do
1179 + jacobi_sub3(phi, col_num, target_i, target_j, i, j + 1, o_row ++ [1], output)
1195 1180 end
1196 1181
1197 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) when (i != target_i or j != target_j) and i != j and j == column_num do
1198 - jacobi_sub3(phi, column_num, target_i, target_j, i + 1, 0, [], output ++ [o_row])
1182 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) when (i != target_i or j != target_j) and i != j and j == col_num do
1183 + jacobi_sub3(phi, col_num, target_i, target_j, i + 1, 0, [], output ++ [o_row])
1199 1184 end
1200 1185
1201 - defp jacobi_sub3(_, column_num, _, _, i, j, _, output) when i == j and j == column_num do
1186 + defp jacobi_sub3(_, col_num, _, _, i, j, _, output) when i == j and j == col_num do
1202 1187 output
1203 1188 end
1204 1189
1205 - defp jacobi_sub3(phi, column_num, target_i, target_j, i, j, o_row, output) do
1206 - jacobi_sub3(phi, column_num, target_i, target_j, i, j + 1, o_row ++ [0], output)
1190 + defp jacobi_sub3(phi, col_num, target_i, target_j, i, j, o_row, output) do
1191 + jacobi_sub3(phi, col_num, target_i, target_j, i, j + 1, o_row ++ [0], output)
1207 1192 end
1208 1193
1209 1194 defp jacobi_sub4({list, index}) do
  @@ -1224,31 +1209,31 @@ defmodule MatrixOperation do
1224 1209 ]
1225 1210 ]
1226 1211 """
1227 - def svd(a, loop_num) do
1212 + def svd(a, iter_num) do
1228 1213 a_t = transpose(a)
1229 - svd_sub(a, a_t, loop_num)
1214 + svd_sub(a, a_t, iter_num)
1230 1215 end
1231 1216
1232 - def svd_sub(a, a_t, loop_num) when length(a) <= length(a_t) do
1217 + def svd_sub(a, a_t, iter_num) when length(a) <= length(a_t) do
1233 1218 # U matrix
1234 1219 aat = product(a, a_t)
1235 - [sv_sq, u] = jacobi(aat, loop_num)
1220 + [sv_sq, u] = jacobi(aat, iter_num)
1236 1221 # V matirx
1237 1222 ata = product(a_t, a)
1238 - [_, v] = jacobi(ata, loop_num)
1223 + [_, v] = jacobi(ata, iter_num)
1239 1224 # Singular value
1240 1225 s = Enum.map(sv_sq, & :math.sqrt(&1))
1241 1226 # A = USV^t
1242 1227 [s, u, v]
1243 1228 end
1244 1229
1245 - def svd_sub(a, a_t, loop_num) do
1230 + def svd_sub(a, a_t, iter_num) do
1246 1231 # U matrix
1247 1232 aat = product(a, a_t)
1248 - [_, u] = jacobi(aat, loop_num)
1233 + [_, u] = jacobi(aat, iter_num)
1249 1234 # V matirx
1250 1235 ata = product(a_t, a)
1251 - [sv_sq, v] = jacobi(ata, loop_num)
1236 + [sv_sq, v] = jacobi(ata, iter_num)
1252 1237 # Singular value
1253 1238 s = Enum.map(sv_sq, & :math.sqrt(&1))
1254 1239 # A = USV^t
  @@ -1263,28 +1248,28 @@ defmodule MatrixOperation do
1263 1248 iex> MatrixOperation.eigenvalue([[6, 1, 1, 1], [1, 7, 1, 1], [1, 1, 8, 1], [1, 1, 1, 9]], 100)
1264 1249 [10.803886359051251, 7.507748705362773, 6.39227529027387, 5.296089645312106]
1265 1250 """
1266 - def eigenvalue(a, loop_num) do
1267 - eigenvalue_sub(a, 0, loop_num)
1251 + def eigenvalue(a, iter_num) do
1252 + eigenvalue_sub(a, 0, iter_num)
1268 1253 end
1269 1254
1270 - defp eigenvalue_sub(a, count, loop_num) when count != loop_num do
1271 - len_a = length(a)
1272 - u = unit_matrix(len_a)
1273 - q_n = qr_for_ev(a, u, len_a, u, 1)
1255 + defp eigenvalue_sub(a, count, iter_num) when count != iter_num do
1256 + matrix_len = length(a)
1257 + u = unit_matrix(matrix_len)
1258 + q_n = qr_for_ev(a, u, matrix_len, u, 1)
1274 1259 a_k = q_n
1275 - |> transpose
1260 + |> transpose()
1276 1261 |> product(a)
1277 1262 |> product(q_n)
1278 - eigenvalue_sub(a_k, count+1, loop_num)
1263 + eigenvalue_sub(a_k, count+1, iter_num)
1279 1264 end
1280 1265
1281 1266 defp eigenvalue_sub(a_k, _, _) do
1282 1267 a_k
1283 - |> Enum.with_index
1268 + |> Enum.with_index()
1284 1269 |> Enum.map(fn {x, i} -> Enum.at(x, i) end)
1285 1270 end
1286 1271
1287 - defp qr_for_ev(a, q, len_a, u, num) when len_a != num do
1272 + defp qr_for_ev(a, q, matrix_len, u, num) when matrix_len != num do
1288 1273 h = get_one_column(a, num)
1289 1274 |> replace_zero(num-1)
1290 1275 |> householder_for_qr(num-1, u)
  @@ -1292,7 +1277,7 @@ defmodule MatrixOperation do
1292 1277 a_n = product(h, a)
1293 1278 q_n = product(q, h)
1294 1279
1295 - qr_for_ev(a_n, q_n, len_a, u, num+1)
1280 + qr_for_ev(a_n, q_n, matrix_len, u, num+1)
1296 1281 end
1297 1282
1298 1283 defp qr_for_ev(_, q_n, _, _, _) do
  @@ -1301,15 +1286,15 @@ defmodule MatrixOperation do
1301 1286
1302 1287 defp replace_zero(list, thresh_num) do
1303 1288 list
1304 - |> Enum.with_index
1289 + |> Enum.with_index()
1305 1290 |> Enum.map(fn {x, i} -> if(i < thresh_num, do: 0, else: x) end)
1306 1291 end
1307 1292
1308 1293 defp householder_for_qr(col, index, u) do
1309 1294 col_norm = col
1310 1295 |> Enum.map(& &1*&1)
1311 - |> Enum.sum
1312 - |> :math.sqrt
1296 + |> Enum.sum()
1297 + |> :math.sqrt()
1313 1298
1314 1299 top = Enum.at(col, index)
1315 1300 top_cn = if(top >= 0, do: top + col_norm, else: top - col_norm)
  @@ -1324,6 +1309,42 @@ defmodule MatrixOperation do
1324 1309 subtract(u, m)
1325 1310 end
1326 1311
1312 + @doc """
1313 + Calculate eigenvector by using QR decomposition
1314 + #### Examples
1315 + iex> MatrixOperation.eigenvector([[1, 4, 5], [4, 2, 6], [5, 6, 3]], 500)
1316 + [
1317 + [0.49659978457191395, 0.577350269219531, 0.6481167492812223],
1318 + [0.31298567717874254, 0.5773502691622936, -0.7541264035190053],
1319 + [-0.8095854617817337, 0.5773502692195656, 0.10600965431255223]
1320 + ]
1321 + """
1322 + def eigenvector(a, iter_num) do
1323 + delta = 0.0001 # avoid division by zero
1324 + a
1325 + |> eigenvalue(iter_num)
1326 + |> Enum.map(
1327 + & eigenvalue_shift(a, -&1+delta)
1328 + |> inverse_matrix()
1329 + |> power_iteration(iter_num)
1330 + |> eigenvector_sub()
1331 + )
1332 + |> Enum.map(& const_multiple(delta, &1))
1333 + end
1334 +
1335 + defp eigenvalue_shift(a, ev) do
1336 + unit = a
1337 + |> length
1338 + |> unit_matrix()
1339 + b = const_multiple(ev, unit)
1340 + add(a, b)
1341 + end
1342 +
1343 + defp eigenvector_sub(a) do
1344 + [_first, second] = a
1345 + second
1346 + end
1347 +
1327 1348 @doc """
1328 1349 Matrix diagonalization
1329 1350 #### Examples
  @@ -1332,9 +1353,9 @@ defmodule MatrixOperation do
1332 1353 iex> MatrixOperation.diagonalization([[2, 1, -1], [1, 1, 0], [-1, 0, 1]], 100)
1333 1354 [[3.000000000000001, 0, 0], [0, 1.0, 0], [0, 0, 0]]
1334 1355 """
1335 - def diagonalization(a, loop_num) do
1336 - eigenvalue(a, loop_num)
1337 - |> diagonalization_condition
1356 + def diagonalization(a, iter_num) do
1357 + eigenvalue(a, iter_num)
1358 + |> diagonalization_condition()
1338 1359 |> Enum.map(& Enum.map(&1, fn x -> zero_approximation(x) end))
1339 1360 end
1340 1361
  @@ -1344,7 +1365,7 @@ defmodule MatrixOperation do
1344 1365
1345 1366 defp diagonalization_condition(a) do
1346 1367 a
1347 - |> Enum.with_index
1368 + |> Enum.with_index()
1348 1369 |> Enum.map(& diagonalization_sub(&1, length(a), 0, []))
1349 1370 end
1350 1371
  @@ -1366,11 +1387,11 @@ defmodule MatrixOperation do
1366 1387 iex> MatrixOperation.singular_value([[1, 2, 3, 1], [2, 4, 1, 5], [3, 3, 10, 8]], 100)
1367 1388 [14.912172620559879, 4.236463407782015, 1.6369134152873956, 0.0]
1368 1389 """
1369 - def singular_value(a, loop_num) do
1390 + def singular_value(a, iter_num) do
1370 1391 a
1371 - |> transpose
1392 + |> transpose()
1372 1393 |> product(a)
1373 - |> eigenvalue(loop_num)
1394 + |> eigenvalue(iter_num)
1374 1395 |> Enum.map(& zero_approximation(&1))
1375 1396 |> Enum.map(& :math.sqrt(&1))
1376 1397 end
  @@ -1383,9 +1404,24 @@ defmodule MatrixOperation do
1383 1404 iex> MatrixOperation.rank([[2, 3, 4, 2], [1, 4, 2, 3], [2, 1, 4, 4]], 100)
1384 1405 3
1385 1406 """
1386 - def rank(matrix, loop_num) do
1387 - singular_value(matrix, loop_num)
1388 - |> count_finite_values
1407 + def rank(matrix, iter_num) do
1408 + matrix
1409 + |> singular_value(iter_num)
1410 + |> count_finite_values()
1411 + end
1412 +
1413 + defp count_finite_values(x) when is_list(x) do
1414 + x
1415 + |> Enum.map(&count_finite_values(&1))
1416 + |> Enum.sum()
1417 + end
1418 +
1419 + defp count_finite_values(x) when is_number(x) and x == 0 do
1420 + 0
1421 + end
1422 +
1423 + defp count_finite_values(x) when is_number(x) do
1424 + 1
1389 1425 end
1390 1426
1391 1427 @doc """
  @@ -1400,8 +1436,8 @@ defmodule MatrixOperation do
1400 1436 a
1401 1437 |> Enum.map(& Enum.map(&1, fn x -> x * x end))
1402 1438 |> Enum.map(& Enum.sum(&1))
1403 - |> Enum.sum
1404 - |> :math.pow(0.5)
1439 + |> Enum.sum()
1440 + |> :math.sqrt()
1405 1441 end
1406 1442
1407 1443 @doc """
  @@ -1416,7 +1452,7 @@ defmodule MatrixOperation do
1416 1452 a
1417 1453 |> Enum.map(& Enum.map(&1, fn x -> if(x > 0, do: x, else: -x) end))
1418 1454 |> Enum.map(& Enum.sum(&1))
1419 - |> Enum.max
1455 + |> Enum.max()
1420 1456 end
1421 1457
1422 1458 @doc """
  @@ -1430,7 +1466,7 @@ defmodule MatrixOperation do
1430 1466 def two_norm(a) do
1431 1467 a
1432 1468 |> singular_value(100)
1433 - |> Enum.max
1469 + |> Enum.max()
1434 1470 end
1435 1471
1436 1472 @doc """
  @@ -1442,10 +1478,11 @@ defmodule MatrixOperation do
1442 1478 10
1443 1479 """
1444 1480 def max_norm(a) do
1445 - transpose(a)
1481 + a
1482 + |> transpose()
1446 1483 |> Enum.map(& Enum.map(&1, fn x -> if(x > 0, do: x, else: -x) end))
1447 1484 |> Enum.map(& Enum.sum(&1))
1448 - |> Enum.max
1485 + |> Enum.max()
1449 1486 end
1450 1487
1451 1488 @doc """
  @@ -1459,7 +1496,7 @@ defmodule MatrixOperation do
1459 1496 """
1460 1497 def variance_covariance_matrix(data) do
1461 1498 x = data
1462 - |> transpose
1499 + |> transpose()
1463 1500 |> Enum.map(& Enum.map(&1, fn x -> x - Enum.sum(&1)/length(&1) end))
1464 1501 xt = transpose(x)
1465 1502 xtx = product(x, xt)
  @@ -4,7 +4,7 @@ defmodule MatrixOperation.MixProject do
4 4 def project do
5 5 [
6 6 app: :matrix_operation,
7 - version: "0.3.5",
7 + version: "0.3.6",
8 8 elixir: "~> 1.7",
9 9 description: "Matrix operation library",
10 10 start_permanent: Mix.env() == :prod,