【问题标题】：Logistic regression in Julia using Optim.jl使用 Optim.jl 在 Julia 中进行逻辑回归
【发布时间】：2015-12-18 14:51:55
【问题描述】：

我正在尝试在 Julia 中实现一个简单的正则化逻辑回归算法。我想使用 Optim.jl 库来最小化我的成本函数，但我无法让它工作。

我的代价函数和梯度如下：

function cost(X, y, theta, lambda)
    m = length(y)
    h = sigmoid(X * theta)
    reg = (lambda / (2*m)) * sum(theta[2:end].^2)
    J = (1/m) * sum( (-y).*log(h) - (1-y).*log(1-h) ) + reg
    return J
end

function grad(X, y, theta, lambda, gradient)
    m = length(y)
    h = sigmoid(X * theta)
    # gradient = zeros(size(theta))
    gradient = (1/m) * X' * (h - y)
    gradient[2:end] = gradient[2:end] + (lambda/m) * theta[2:end]
    return gradient
end

（其中theta 是假设函数的参数向量，lambda 是正则化参数。）

然后，根据这里给出的说明：https://github.com/JuliaOpt/Optim.jl我尝试这样调用优化函数：

# those are handle functions I define to pass them as arguments:
c(theta::Vector) = cost(X, y, theta, lambda)
g!(theta::Vector, gradient::Vector) = grad(X, y, theta, lambda, gradient)

# then I do
optimize(c,some_initial_theta) 
# or maybe
optimize(c,g!,initial_theta,method = :l_bfgs) # try a different algorithm

在这两种情况下，它都表示它无法收敛，并且输出看起来有点尴尬：

julia> optimize(c,initial_theta)
Results of Optimization Algorithm
 * Algorithm: Nelder-Mead
 * Starting Point: [0.0,0.0,0.0,0.0,0.0]
 * Minimum: [1.7787162051775145,3.4584135105727145,-6.659680628594007,4.776952006060713,1.5034743945407143]
 * Value of Function at Minimum: -Inf
 * Iterations: 1000
 * Convergence: false
   * |x - x'| < NaN: false
   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false
   * |g(x)| < NaN: false
   * Exceeded Maximum Number of Iterations: true
 * Objective Function Calls: 1013
 * Gradient Call: 0

julia> optimize(c,g!,initial_theta,method = :l_bfgs)
Results of Optimization Algorithm
 * Algorithm: L-BFGS
 * Starting Point: [0.0,0.0,0.0,0.0,0.0]
 * Minimum: [-6.7055e-320,-2.235e-320,-6.7055e-320,-2.244e-320,-6.339759952602652e-7]
 * Value of Function at Minimum: 0.693148
 * Iterations: 1
 * Convergence: false
   * |x - x'| < 1.0e-32: false
   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false
   * |g(x)| < 1.0e-08: false
   * Exceeded Maximum Number of Iterations: false
 * Objective Function Calls: 75
 * Gradient Call: 75

问题

我的方法（来自我的第一个代码清单）不正确吗？还是我滥用了 Optim.jl 函数？无论哪种方式，在这里定义和最小化成本函数的正确方法是什么？

这是我第一次接触 Julia，可能我做错了什么，但我不知道具体是什么。任何帮助将不胜感激！

编辑

X 和 y 是训练集，X 是一个 90x5 矩阵，y 一个 90x1 向量（也就是说，我的训练集取自 Iris - 我不知道不认为这很重要）。

【问题讨论】：

编程问题在这里是题外话。在 stack-overflow 上问这个问题你可能会有更好的运气。
感谢@MarcClaesen 的回复，我在这里问它是因为我不确定我的方法是否正确，而不仅仅是实现。无论如何，我会在 SO 尝试。
如果您的问题中隐藏了方法问题，您可能需要对其进行编辑以将其放在首位。然后，编程问题就不那么突出了，可以在此处考虑主题。但是，如果您想尝试 SO，请将帖子标记为移动而不是重新发布。谢谢！
我尽可能多地澄清了我的问题，但也将其标记为转移到 SO。感谢您的建议。
我没有看到任何控制输出级别的选项，除了 show_trace，它可能会显示有用的信息并帮助您了解正在发生的事情。也许尝试不同的起始值会有所帮助。默认公差是疯狂的（不是一个好兆头）。但这是关键，您使用的包并不完全是最重要的。每github.com/JuliaOpt/Optim.jl ，“因为它是从头开始开发的，所以它不如基于 C 的 NLOpt 包强大。对于准确性必须不容置疑的工作，我们建议使用 NLOpt 包。”使用更好的东西，无论是在 Julia 还是在其他地方。

标签： regression julia mathematical-optimization

【解决方案1】：

这是一个使用 Optim.jl 的自分化功能的非正则逻辑回归示例。它可能对您自己的实现有所帮助。

using Optim

const X = rand(100, 3)
const true_β = [5,2,4]
const true_y =  1 ./ (1 + exp(-X*true_β))

function objective(β)
    y = 1 ./ (1 + exp(-X*β))
    return sum((y - true_y).^2)  # Use SSE, non-standard for log. reg.
end

println(optimize(objective, [3.0,3.0,3.0],
                autodiff=true, method=LBFGS()))

这给了我

Results of Optimization Algorithm
 * Algorithm: L-BFGS
 * Starting Point: [3.0,3.0,3.0]
 * Minimizer: [4.999999945789497,1.9999999853962256,4.0000000047769495]
 * Minimum: 0.000000
 * Iterations: 14
 * Convergence: true
   * |x - x'| < 1.0e-32: false
   * |f(x) - f(x')| / |f(x)| < 1.0e-08: false
   * |g(x)| < 1.0e-08: true
   * Exceeded Maximum Number of Iterations: false
 * Objective Function Calls: 53
 * Gradient Call: 53

【讨论】：

使用这种方法有什么缺点吗？使用自动微分不是比显式提供梯度时慢吗？
它较慢（通常），但可能有助于检查您的实施（我怀疑这可能是您遇到问题的原因）。
你能更新一下这个脚本吗？不适用于 Julia v1.2。

【解决方案2】：

您可以在下面找到我使用闭包和柯里化进行逻辑回归的成本和梯度计算函数（适用于那些习惯于返回成本和梯度的函数的人的版本）：

function cost_gradient(θ, X, y, λ)
    m = length(y)
    return (θ::Array) -> begin 
        h = sigmoid(X * θ)   
        J = (1 / m) * sum(-y .* log(h) .- (1 - y) .* log(1 - h)) + λ / (2 * m) * sum(θ[2:end] .^ 2)         
    end, (θ::Array, storage::Array) -> begin  
        h = sigmoid(X * θ) 
        storage[:] = (1 / m) * (X' * (h .- y)) + (λ / m) * [0; θ[2:end]]        
    end
end

Sigmoid函数实现：

sigmoid(z) = 1.0 ./ (1.0 + exp(-z))

要在 Optim.jl 中应用 cost_gradient，请执行以下操作：

using Optim
#...
# Prerequisites:
# X size is (m,d), where d is the number of training set features
# y size is (m,1)
# λ as the regularization parameter, e.g 1.5
# ITERATIONS number of iterations, e.g. 1000
X=[ones(size(X,1)) X] #add x_0=1.0 column; now X size is (m,d+1)
initialθ = zeros(size(X,2),1) #initialTheta size is (d+1, 1)
cost, gradient! = cost_gradient(initialθ, X, y, λ)
res = optimize(cost, gradient!, initialθ, method = ConjugateGradient(), iterations = ITERATIONS);
θ = Optim.minimizer(res);

现在，您可以轻松预测（例如训练集验证）：

predictions = sigmoid(X * θ) #X size is (m,d+1)

要么尝试我的方法，要么将其与您的实现进行比较。

【讨论】：