Merge e51f9d4 into aa93a76

lib/matrex/algorithms.ex

@@ -272,22 +272,24 @@ defmodule Matrex.Algorithms do
 
     cond do
       data.f2 in [:nan, :inf, :neg_inf] ->
+        # IO.puts("iter2: inf ")
         {false, data}
 
       data.f2 > data.f1 + data.z1 * @rho * data.d1 or data.d2 > -@sig * data.d1 ->
+        # IO.puts("iter2: first ")
         {false, data}
 
       data.d2 > @sig * data.d1 ->
-        # IO.puts("second #{data.d1}, #{data.d2}")
+        # IO.puts("iter2: second #{data.d1}, #{data.d2}")
         {true, data}
 
       data.m == 0 ->
-        # IO.puts("third")
+        # IO.puts("iter2: third")
         {false, data}
 
       true ->
         z2 = z2(data, data.limit)
-        # IO.puts("none, z2=#{z2}")
+        # IO.puts("iter2: none, z2=#{z2}")
 
         data = %{
           data
@@ -401,12 +403,12 @@ defmodule Matrex.Algorithms do
   end
 
   @doc """
-  Linear regression cost and gradient function with regularization from Andrew Ng's course (ex3).
+  Logistic regression cost and gradient function with regularization from Andrew Ng's course (ex3).
 
   Computes the cost of using `theta` as the parameter for regularized logistic regression and the
   gradient of the cost w.r.t. to the parameters.
 
-  Compatible with `fmincg/4` algorithm from thise module.
+  Compatible with `fmincg/4` algorithm from this module.
 
   `theta`  — parameters, to compute cost for
 
@@ -583,6 +585,99 @@ defmodule Matrex.Algorithms do
     accuracy
   end
 
+  @doc """
+  Linear regression cost and gradient function, with no normalization.
+
+  Computes the cost of using `theta` as the parameter for regularized linear regression and the
+  gradient of the cost w.r.t. to the parameters.
+
+  Compatible with `fmincg/4` algorithm from thise module.
+
+  `theta`  — parameters, to compute cost for
+
+  `X`  — training data input.
+
+  `y`  — training data output.
+
+  `lambda`  — regularization parameter.
+
+  """
+  @spec linear_cost_fun(Matrex.t(), {Matrex.t(), Matrex.t(), number, non_neg_integer}, pos_integer) ::
+          {float, Matrex.t()}
+  def linear_cost_fun(
+        %Matrex{} = theta,
+        {%Matrex{} = x, %Matrex{} = y, lambda} = _params,
+        _iteration \\ 0
+      ) when is_number(lambda) do
+    n = Enum.count(y)
+    h! = x |> Matrex.dot(theta)
+
+    h_sub_y = Matrex.subtract(h!,y)
+
+    j = 1.0/(2*n) *
+      (Matrex.dot( h_sub_y |> Matrex.transpose, h_sub_y) |> Matrex.scalar())
+
+    grad = x |> Matrex.transpose |> Matrex.dot(h_sub_y) |> Matrex.apply(& &1 * lambda / n)
+
+    {j, grad}
+  end
+
+  @doc """
+  Fit polynomial function based on given `x` and `y` using gradient descent (`fmincg/4`) using
+  `linear_cost_fun/4`.
+
+  This approach produces decent results for many datasets. It isn't as efficient as using
+  least squares, but is still useful and provides a goode example of how to optimize general
+  functions.
+
+  Note that gradient descent won't always converge well for polynomials with linear cost function.
+  If this happens for your dataset try adjusting `opts` parameters. Uses the `fmincg/4` algorithm
+  from this module.
+
+  `x`  — training data input.
+
+  `y`  — training data output.
+
+  `opts` - algorithm parameters
+    `lambda`  — regularization parameter.
+    `iterations`  — regularization parameter.
+
+  """
+  @spec fit_poly(Matrex.t(), Matrex.t(), pos_integer, keyword() ) ::
+      %{ coefs: keyword(), error: float(), fun: (Matrex.t() -> Matrex.t()) }
+  def fit_poly(x, y, degree, opts \\ []) do
+    iterations = Keyword.get(opts, :iterations, 100)
+    lambda  = Keyword.get(opts, :lambda, 1.0)
+
+    {m, n} = Matrex.size(y)
+    unless m >= n, do: raise %ArgumentError{message: "y shape (m,n) must have m > n"}
+
+    xx =
+      for i <- 0..degree, into: [] do
+        x |> Matrex.apply(&:math.pow(&1, i))
+      end |> Matrex.concat()
+
+    theta = Matrex.zeros(degree + 1, 1)
+
+    {sX, fX, _i} =  fmincg(&linear_cost_fun/3, theta, {xx, y, lambda}, iterations)
+
+    coefs = sX |> Enum.to_list() |> Enum.with_index(0) |> Enum.map(fn {x,y} -> {y,x} end)
+    %{coefs: coefs, fun: &poly_func(&1, coefs), error: fX |> Enum.at(-1)}
+  end
+
+  def fit_linear(x, y, opts \\ []) do
+    fit_poly(x, y, 1, opts)
+  end
+
+  defp poly_func(x, coefs) when is_list(coefs) do
+    # coefs_idx = Enum.with_index(coefs, 0)
+    x |> Enum.map(fn x ->
+      Enum.reduce(coefs, 0.0, fn {i, f}, acc ->
+        acc + f * :math.pow(x, i)
+      end)
+    end)
+  end
+
   @doc """
   Function of a surface with two hills.
   """

test/algorithms_test.exs

@@ -115,4 +115,57 @@ defmodule AlgorithmsTest do
 
     {x_train, y_train, x_test, y_test}
   end
+
+
+  test "#linear_cost_fun computes cost" do
+    m = Matrex.load("test/rand_array.mtx")
+    y_t = m |> Matrex.submatrix(1..41, 2..2)
+
+    # for linear func, must add `ones` for the offset constant
+    x = m |> Matrex.submatrix(1..41, 1..1)
+    x_t = Matrex.concat(Matrex.ones(Matrex.size(x)), x)
+
+    lambda_t = 0.01
+    theta_t = Matrex.zeros(2, 1)
+
+    expected_j = 5238.50381097561
+
+    expected_grad =
+      Matrex.new(
+        " -0.91246 ; -2.41489 "
+      )
+
+    {j, grad} = Algorithms.linear_cost_fun(theta_t, {x_t, y_t, lambda_t})
+
+    assert grad |> Matrex.subtract(expected_grad) |> Matrex.sum() < 5.0e-6
+    assert j == expected_j
+
+  end
+
+  test "#fit_poly " do
+    m = Matrex.load("test/rand_array.mtx")
+    y = m |> Matrex.submatrix(1..41, 2..2)
+    x = m |> Matrex.submatrix(1..41, 1..1)
+
+    fit = Algorithms.fit_poly(x, y, 2)
+
+    expected_fit = %{
+      coefs: [
+        {0, 37.48050308227539},
+        {1, 6.260676383972168},
+        {2, 6.991103172302246}
+      ],
+      error: 149.0388957698171,
+    }
+
+    # IO.inspect(fit, label: :fit)
+    expected_coefs = expected_fit[:coefs] |> coefs_nums()
+    coefs = fit[:coefs] |> coefs_nums()
+
+    assert coefs |> Matrex.subtract(expected_coefs) |> Matrex.sum() < 1.0e-5
+  end
+
+  defp coefs_nums(c) do
+    [c |> Enum.map(& &1 |> elem(1))] |> Matrex.new()
+  end
 end