In [14]:
using Plots, LinearAlgebra, Random, SparseArrays
using Optim
In [15]:
# example of Rosenbrock function 
# NOTE, This doesn't work, see the next cell! 
# See these things are hard to use :) 
function f(x) 
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
function g!(x::Vector, storage::Vector)
storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
storage[2] = 200.0 * (x[2] - x[1]^2)
end

soln = optimize(f, g!, [0.0, 0.0], GradientDescent())
Out[15]:
 * Status: failure

 * Candidate solution
    Final objective value:     NaN

 * Found with
    Algorithm:     Gradient Descent

 * Convergence measures
    |x - x'|               = NaN ≰ 0.0e+00
    |x - x'|/|x'|          = NaN ≰ 0.0e+00
    |f(x) - f(x')|         = NaN ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = NaN ≰ 0.0e+00
    |g(x)|                 = NaN ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    0
    f(x) calls:    1
    ∇f(x) calls:   1
In [16]:
# example of Rosenbrock function
function f(x) 
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
function g!(storage::Vector, x::Vector)
storage[1] = -2.0 .* (1.0 - x[1]) - 400.0 .* (x[2] - x[1]^2) * x[1]
storage[2] = 200.0 .* (x[2] - x[1]^2)
end

soln = optimize(f, g!, [0.0, 0.0], BFGS())
Out[16]:
 * Status: success

 * Candidate solution
    Final objective value:     7.645688e-21

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 3.48e-07 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.48e-07 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.91e-14 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 9.03e+06 ≰ 0.0e+00
    |g(x)|                 = 2.32e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    16
    f(x) calls:    53
    ∇f(x) calls:   53
In [ ]:
In [9]:
soln = optimize(f, g!, [0.0, 0.0], GradientDescent())
Out[9]:
 * Status: failure (reached maximum number of iterations)

 * Candidate solution
    Final objective value:     4.154719e-03

 * Found with
    Algorithm:     Gradient Descent

 * Convergence measures
    |x - x'|               = 1.82e-04 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.95e-04 ≰ 0.0e+00
    |f(x) - f(x')|         = 8.18e-06 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.97e-03 ≰ 0.0e+00
    |g(x)|                 = 8.21e-02 ≰ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    1000
    f(x) calls:    2532
    ∇f(x) calls:   2532
In [10]:
# now, we do the matrix factorization example
# originally from Poblano example2

function matrix_approx_function(x::Vector, A::Matrix, r::Int)
    # unpack U and V from x
    m,n = size(A)
    U = reshape(x[1:m*r],m,r)
    V = reshape(x[(m*r+1):end],n,r)
    return 0.5*norm(A - U*V')^2 
end

function matrix_approx_gradient!(storage::Vector, x::Vector, A::Matrix, r::Int)
    m,n = size(A)
    U = reshape(x[1:m*r],m,r)
    V = reshape(x[(m*r+1):end],n,r)
    D = A - U*V'
    storage[1:(m*r)] = -vec(D*V)
    storage[(m*r+1):end] = -vec(D'*U)
end
Out[10]:
matrix_approx_gradient! (generic function with 1 method)
In [19]:
m = 10
n = 8
#A = randn(m,n)
A = Matrix(1.0I,m,n)
r = 2
myf = x -> matrix_approx_function(x, A, r)
myg! = (x, storage) -> matrix_approx_gradient!(x, storage, A, r)

soln = optimize(myf, myg!, ones(m*r+n*r), BFGS(), Optim.Options(f_tol = 1e-8))
#soln = optimize(myf, myg!, randn(m*r+n*r), BFGS(), Optim.Options(f_tol = 1e-8))
x = Optim.minimizer(soln)
@show soln
Uopt = reshape(x[(1:m*r)],m,r)
Vopt = reshape(x[(m*r+1):end],n,r)
objval = 2*myf(x)
opterr = norm(A-Uopt*Vopt')^2

Usvd,Ssvd,Vsvd = svd(A)
svderr = norm(A-Usvd[:,1:r]*Diagonal(Ssvd[1:r])*Vsvd[:,1:r]')^2
@show objval
@show opterr
@show svderr
; # hide final output in JuliaBox

soln =  * Status: success

 * Candidate solution
    Final objective value:     3.500000e+00

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 4.60e-06 ≰ 0.0e+00
    |x - x'|/|x'|          = 3.64e-06 ≰ 0.0e+00
    |f(x) - f(x')|         = 4.22e-11 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.21e-11 ≤ 1.0e-08
    |g(x)|                 = 5.43e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    11
    f(x) calls:    31
    ∇f(x) calls:   31

objval = 7.000000000000001
opterr = 7.000000000000001
svderr = 5.999999999999999
In [20]:
Uopt
Out[20]:
10×2 Matrix{Float64}:
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
  1.26522     1.26522
 -5.0131e-8  -5.0131e-8
 -5.0131e-8  -5.0131e-8
In [13]:
Vopt
Out[13]:
8×2 Matrix{Float64}:
  0.0659192  -0.435066
 -0.341172   -0.270793
 -0.311504    0.552242
  0.173043   -0.596362
 -0.225341   -0.256013
 -0.046164   -0.396231
 -0.269815   -0.463872
  0.817992    0.289806
In [ ]: