JuMP write_to_file: Writing legacy model to file is much smaller

Specific Domains Optimization (Mathematical)
1

I construct a model for a nonlinear optimization problem in the Julia JuMP code below, using either the non-legacy (use_legacy_model = false) or legacy (use_legacy_model = true) model interface and save the model to a .nl file using write_to_file. The .nl file created from the legacy model interface is much smaller (e.g. 229KB) than the .nl created from the non-legacy model interface (e.g. 4.5MB). What is the reason for the significantly different file sizes?

-rw-rw-r-- 1 stuart stuart 229K Jun 6 13:43 model_legacy.nl
-rw-rw-r-- 1 stuart stuart 4.5M Jun 6 13:44 model_non_legacy.nl

using JuMP

use_legacy_model = false # Whether to use legacy JuMP nonlinear constraints, expressions, and objectives.

N1 = 64
N = 64
zpad = 2
M = zpad*N

s = rand(M,1) 
Br = rand(M, N1)
Bi = rand(M, N1)

Mt = trunc(Int,floor(M/2))+1
Br = Br[1:Mt,:]
Bi = Bi[1:Mt,:]
s = s[1:Mt]
Me = Mt

# Construct the JuMP model.
model_construction_time = @elapsed begin # Start elapsed.

model = Model()

@variable(model, y[1:N1], Bin)

if(~use_legacy_model)
        print("Constructing non-legacy model.\n")

        @expression(model, yh, y.-.5)
        #@expression(model, yh[n=1:N1], y[n]-.5)

        @expression(model, spec_r, Br*yh)
        @expression(model, spec_i, Bi*yh)

        #@expression(model, spec_sq_mag, spec_r.^2 .+ spec_i.^2) # Slow
        @expression(model, spec_sq_mag[m=1:Me], spec_r[m]^2 + spec_i[m]^2) # Fast

        #@expression(model, gamma, sum((spec_sq_mag.-s).^2)) # Slow
        @expression(model, gamma, sum((spec_sq_mag[m]-s[m])^2 for m in 1:Me)) # Fast

        @objective(model, Min, gamma)
else
        print("Constructing legacy model.\n")

        @expression(model, yh[n=1:N1], y[n]-.5)
        @expression(model, spec_r, Br*yh)
        @expression(model, spec_i, Bi*yh)

        @NLexpression(model, spec_sq_mag[m=1:Me], spec_r[m]^2 + spec_i[m]^2)

        @variable(model, gamma >= 0)

        @NLconstraint(model, gamma >= sum((spec_sq_mag[m]-s[m])^2 for m in 1:Me))
        @objective(model, Min, gamma)
end

end # End elapsed.
println("\nJuMP model construction time: $(model_construction_time) [s].")

# Write the JuMP model to a .nl file.
model_write_time = @elapsed begin # Start elapsed.

write_to_file(model,"model.nl")

end # End elapsed.
println("\nJuMP model write time: $(model_write_time) [s].")
2

The issue is quite subtle, and is because of:

@expression(model, spec_sq_mag[m=1:Me], spec_r[m]^2 + spec_i[m]^2)

This line will automatically try to create QuadExpr objects.

In most cases, this is the right thing to do, because it allows solvers to specialize on quadratic constraints and objectives. However, when writing to a .nl file, this has a bad outcome.

Consider this simple example:

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> y = sum(x)
x[1] + x[2]

julia> y^2
x[1]² + 2 x[2]*x[1] + x[2]²

julia> @force_nonlinear(y^2)
(x[1] + x[2]) ^ 2

The .nl file writes things out as a tape of the expression graph.

If we write out the QuadExpr y^2 as a tape, it requires 11 nodes:

[+, *, x[1], x[2], *, 2, x[2], x[1], *, x[2], x[2]]

If we write out the NonlinearExpr, it requires only 5 nodes:

[^, +, x[1], x[2], 2]

You can fix things for your code with:

@expression(
    model, 
    spec_sq_mag[m=1:Me],
    @force_nonlinear(spec_r[m]^2 + spec_i[m]^2),
)

See: JuMP · JuMP

3

If @force_nonlinear is used, will the solver still recognize the constraint or expression as being quadratic? If not, then it seems that @force_nonlinear should only be used if the model must be written to a file and the file size must be small.

4

will the solver still recognize the constraint or expression as being quadratic

No.

it seems that @force_nonlinear should only be used if the model must be written to a file and the file size must be small

Yes. This is why @force_nonlinear is not the default, and it is one reason that we added it.

You can do better by slightly reformulating your problem:

using JuMP
N1, N, zpad = 64, 64, 2
M = zpad * N
Mt = ceil(Int, M / 2)
s, Br, Bi = rand(Mt), rand(Mt, N1), rand(Mt, N1);
model = Model();
@variable(model, y[1:N1], Bin);
@expression(model, yh, y .- 0.5);
@variable(model, spec_r[1:Mt])
@variable(model, spec_i[1:Mt])
@expression(model, spec_r .== Br * yh)
@expression(model, spec_i .== Bi * yh)
@expression(model, spec_sq_mag[m in 1:Mt], spec_r[m]^2 + spec_i[m]^2)
@expression(model, gamma, sum((spec_sq_mag[m]-s[m])^2 for m in 1:Mt))
@objective(model, Min, gamma)
write_to_file(model, "model.nl")

Now it is only 207K:

-rw-r--r--  1 oscar  staff   207K Jun  7 12:06 model.nl
5

Thanks for your suggestions. Your formula for Mt is not the same as in my original code. ceil(Int, M / 2) does not equal floor(Int, M / 2)+1 when Mt is even.

1 Like
6

Indeed :grimacing: I guess I didn’t think carefully enough

7

When I ran your JuMP code, the resulting .nl file is only 6.7 KB, rather than 207 KB!

-rw-rw-r-- 1 stuart stuart 6.7K Jun  6 17:48 model.nl
8

I copied the wrong thing:

using JuMP
N1, N, zpad = 64, 64, 2
M = zpad * N
Mt = ceil(Int, M / 2)
s, Br, Bi = rand(Mt), rand(Mt, N1), rand(Mt, N1);
model = Model();
@variable(model, y[1:N1], Bin);
@expression(model, yh, y .- 0.5);
@variable(model, spec_r[1:Mt])
@variable(model, spec_i[1:Mt])
@constraint(model, spec_r .== Br * yh)
@constraint(model, spec_i .== Bi * yh)
@expression(model, spec_sq_mag[m in 1:Mt], spec_r[m]^2 + spec_i[m]^2)
@expression(model, gamma, sum((spec_sq_mag[m]-s[m])^2 for m in 1:Mt))
@objective(model, Min, gamma)
write_to_file(model, "model.nl")

Note @constraint(model, spec_r .== Br * yh)

1 Like