New to Optimization - Guidance

Hi,

I’m trying to implement an algorithm from a paper[1] which includes a nonlinear optimization. I’ve never done any optimization before, aside from fitting curves with inverting matrices. I’m struggling to work out what goes where and what I’m allowed to do with things in JuMP.

The quantity to minimise is a sum of squares of a function of some input data. The function involves several nested nonlinear operations on the optimisation variables. I tried implementing the whole thing using a closure to enable access to the input data. This results in an error in a macro parsing step:

ERROR: LoadError: MethodError: no method matching parseNLExpr_runtime(::JuMP.Mod el, ::Array{JuMP.Variable,2}, ::Array{ReverseDiffSparse.NodeData,1}, ::Int64, :: Array{Float64,1})

where the second variable is not matched. The best I could guess from this is that the objective function can only be a function of a single optimization variable, but this seems like an unlikely limitation.

At one point which I put all the variables in one array I got this error:

ERROR: LoadError: Unexpected vector v[i] free ∀ i ∈ {1,2,…,23,24} in nonlinear expression. Nonlinear expressions may contain only scalar expressions.

In this case, the vector of variables v was always ‘dereferenced’ with a getindex. Is this not allowed?

I don’t understand from reading the documentation what kinds of function should go in constraints, objectives, and expressions. My questions:

Should I expect to be able to write a closure to pass to @NLobjective which is a function of 20ish variables, and contains an array of 40 or so floats which are the input data? Or do I write lots of the sub-computations as @constraints setting intermediate variables, and then writing an @NLobjective which just squares and sums an intermediate array? Are arrays allowed in Nonlinear solves?

[1] http://www.robots.ox.ac.uk/~dclaus/publications/claus05plumbline.pdf

Can you post a short code snippet that reproduces the errors? It’s really
difficult to diagnose these problems otherwise.

-Joey

I’m afraid I have no idea how to reduce this to a snippet. I’ll post a script which reproduces the error, for anyone kind enough to have a look. Hmm, no code files as attachments. Inline it is then. Apologies for the huge dump. The function findlenscalibration and the closure generaterationallenscostfn contain the things I’m trying to understand.

using JuMP, Parameters

import Base: dot

@with_kw immutable LensCalibSolveParams
    a_init = 1.0
    b_init = 2*10^7
    w_init = 0.5
    imagewidth = 1920
    imageheight = 1080
    maxedgelsperline = 20

end

@with_kw immutable Conic
    a = 0.0
    b = 0.0
    c = 0.0
    d = 0.0
    e = 0.0
    f = 0.0
end

function dot{T<:AbstractArray}(C::Conic, b::T)
    return C.a*b[1]+C.b*b[2]+C.c*b[3]+C.d*b[4]+C.e*b[5]+C.f*b[6]
end

function conicdesignmatrix!(matrix, xs, ys)
    for r in 1:size(matrix,1)
        x = xs[r]
        y = ys[r]
        matrix[r,1] = x^2
        matrix[r,2] = x*y
        matrix[r,3] = y^2
        matrix[r,4] = x
        matrix[r,5] = y
        matrix[r,6] = 1.0
    end
    matrix
end

function fitconic(xs, ys, design_matrix)
    conicdesignmatrix!(design_matrix, xs, ys)
    U, S, Vs = svd(design_matrix, thin=false)
    # The null space is the column corresponding to the smallest singular value.
    # The SVD returns the singular values sorted, so this is just the final column of V:
    Conic(Vs'[:,6]...)
end



"A closure to create the cost function"
function generaterationallenscostfn(params, linesxs, linesys)
    @unpack_LensCalibSolveParams params

    # allocate memory
    design_matrix = zeros(5,6)

    # Fit a conic to each line, then find the dot product of the conic
    # with each corresponding lifted point. Store the dot products in an array.
    # This needs a 2D array with dimensions NoOfLinesxMaxNoOfEdgelsPerLine
    conics = Array(Conic, length(linesxs))
    for i in 1:length(linesxs)
        conics[i] = fitconic(linesxs[i], linesys[i], design_matrix)
    end

    theta_dot_chis = zeros(length(linesxs), maxedgelsperline)
    for j in 1:length(linesxs), i in 1:min(maxedgelsperline, length(linesxs[j]))
        x = linesxs[j][i]
        y = linesys[j][i]
        chi = [x^2, x*y, y^2, x, y, 1]
        theta_dot_chis[j,i] = dot(conics[j], chi)
    end

    nlines = length(linesxs)

    A = zeros(3, 6)
    Axx = 0.0
    Axy = 0.0
    Ax = 0.0
    Ayy = 0.0
    Ay = 0.0

    return function(v) # v = [a, b, w, lines...]
        result = 0.0

        asquared = v[1]^2
        wsquared = v[3]^2
        A31 = 4*asquared*wsquared

        # Each iteration:
        # For each line:
        #for (lidx, line) in enumerate(lines)
        for linen in 1:nlines
            line = v[linen*3+1:linen*3+3] # 4, 7, 10, 13...
            Axx = sum(line[1]*1.0, line[2]*1.0, line[3]*A31) #dot(A[:,1], line)
            Axy = 0.0#dot(A[:,2], line)
            Ayy = sum(line[1]/asquared, line[2]/asquared, line[3]*4*wsquared)#dot(A[:,3], line)
            Ax = sum((-2*v[2]-imagewidth)*line[1], -imagewidth*line[2], -A31*imagewidth*line[3])#dot(A[:,4], line)
            Ay = sum(line[2]*(-2*v[2]/asquared-imageheight/asquared), line[3]*(-4*wsquared*imageheight))#dot(A[:,5], line)
        # Find A_xx, A_xy... for the line
        # For each point on the line:
            for pidx in 1:nlines
        # Take the conic-point product for this point, and divide it by the denom made
        # up using this line and A_xx A_xy...
        # Accumulate these & return.
                x = linesxs[lidx][pidx]
                y = linesys[lidx][pidx]
                denom = sqrt((2Axx*x+Axy*y+Ax)^2+(2Ayy*y+Axy*x+Ay)^2)
                result += (theta_dot_chis[lidx, pidx]/denom)^2
            end
        end
        return result
    end

end

function initline(xs, ys)
    # Fit a straight line to the points (x,y) in xs, ys to initialise.
    A = zeros(2,2)
    b = zeros(2)
    c = zeros(2)
    A[1,1] = length(xs)
    # BLASification
    A[1,2] = sum(xs)
    A[2,2] = BLAS.dot(length(xs), xs, 1, xs, 1)
    A[2,1] = A[1,2]
    b[1] = sum(ys)
    b[2] = BLAS.dot(length(xs), xs, 1, ys, 1)
    try
        c = inv(A)*b
    catch
        c = pinv(A)*b
    end
    [c... 1]

end


# Computes the initial estimates of straight lines for each of the curves.
function init_lines(linesxs, linesys, params)
    @unpack_LensCalibSolveParams params
    lines = zeros(3,length(linesxs))
    for i in 1:length(linesxs)
        lines[:,i] = initline(linesxs[i], linesys[i])
    end
    lines
end


function findlenscalibration(linesxs, linesys, params)

    @unpack_LensCalibSolveParams params

    nlines = length(linesxs)

    m = Model()

    # estimate initial values:
    lines_init = init_lines(linesxs, linesys, params)
    lines_init = reshape(lines_init, 3*nlines)
    initvals = [a_init b_init w_init lines_init...]

    @variable(m, v[i=1:length(initvals)], start=initvals[i])

    # Get a function to compute the cost, given the constant parameters.
    rationallenscost = generaterationallenscostfn(params, linesxs, linesys)

    JuMP.register(m, :rationallenscost, 1, rationallenscost, autodiff=true)

    @NLobjective(m, Min, rationallenscost(v))

    status = solve(m)

    computeA(getvalue(a), getvalue(b), getvalue(w))
end

line1x = [403.5, 641.5, 926.5, 1220.0, 1476.5]
line1y = [201.0, 169.5, 152.0, 156.5, 179.5]

line2x = [391.0, 633.0, 928.0, 1232.5, 1495.0]
line2y = [465.0, 456.0, 450.0, 449.5, 454.0]

line3x = [405.0, 643.5, 931.0, 1228.5, 1486.0]
line3y = [733.5, 749.5, 756.5, 750.5, 734.5]

# "Vertical" lines:
line4x = [515.0, 507.0, 504.0, 508.0, 516.5, 531.0]
line4y = [184.0, 319.0, 460.0, 603.0, 742.0, 873.0]

line5x  =[779.5, 776.0, 776.0, 778.5, 783.0, 790.0]
line5y = [158.0, 301.0, 452.5, 607.0, 754.5, 893.0]

line6x = [1074.0, 1079.0, 1081.5, 1083.0, 1081.5, 1078.5]
line6y = [151.0, 295.0, 449.0, 605.0, 755.5, 895.5]

line7x = [1355.0, 1365.5, 1371.5, 1370.5, 1364.5, 1352.5]
line7y = [166.0, 305.0, 451.0, 600.0, 743.5, 879.5]

# Collect these together into arrays of arrays:
linesxs = [line1x, line2x, line3x, line4x, line5x, line6x, line7x]
linesys = [line1y, line2y, line3y, line4y, line5y, line6y, line7y]

params = LensCalibSolveParams()

findlenscalibration(linesxs, linesys, params)

As the second error message suggests, only scalar-valued user defined functions are allowed. Instead of passing in a vector of variables v, you’ll need to rewrite the logic so that each argument is a single decision variable (e.g. v[1]).

Unconstrained optimization of a single user-defined function is not a good use case for JuMP; JuMP does nothing for you in this case. You should look at Optim instead.

Thanks Miles. As someone new to the whole domain of optimization, that wasn’t obvious. A question: The lines in v should in fact be normalised to the range -1…1. Is this sufficient constraint to gain benefit from JuMP, or should I still investigate Optim?