google-research

1
import FirstOrderLp, Plots, CSV
2
using LaTeXStrings, DataFrames, Test
3
include("../src/PrimalDualMethods.jl")
4

5
simple_lp1 = FirstOrderLp.linear_programming_problem(
6
  [-Inf],  # variable_lower_bound
7
  [Inf],  # variable_upper_bound
8
  [0.0],  # objective_vector
9
  0.0,                 # objective_constant
10
  ones(1, 1),           # constraint_matrix
11
  [0.0],      # right_hand_side
12
  1,                     # num_equalities
13
)
14
simple_lp2 = FirstOrderLp.linear_programming_problem(
15
  [0.0, 0.0],  # variable_lower_bound
16
  [Inf, Inf],  # variable_upper_bound
17
  [2.0, 1.0],  # objective_vector
18
  0.0,                 # objective_constant
19
  ones(1, 2),           # constraint_matrix
20
  [1.0],      # right_hand_side
21
  1,                     # num_equalities
22
)
23

24
A = [
25
  -1.03324 -0.449087 0.322969 1.36308 1.46373 -1.63236
26
  -0.326885 0.471807 1.73327 -0.505794 0.543352 -1.2965
27
]
28
lp1_primal_optimal_solution = [0.0; 0.0; 0.0; 0.0; 1.0; 1.0]
29
b = A * lp1_primal_optimal_solution
30
lp1 = FirstOrderLp.linear_programming_problem(
31
  zeros(6),  # variable_lower_bound
32
  Inf * ones(6),  # variable_upper_bound
33
  [10.0, 10.0, 10.0, 10.0, 1.0, 1.0],  # objective_vector
34
  0.0,                 # objective_constant
35
  A,           # constraint_matrix
36
  b,      # right_hand_side
37
  2,      # num_equalities
38
)
39

40
A_B = A[:, end-1:end]
41
lp1_dual_optimal_solution = A_B' \ lp1.objective_vector[end-1:end]
42

43
restart_length = 25
44

45
for method_name in ["PDHG", "extragradient", "ADMM"]
46
  @testset "$method_name" begin
47
    for primal_weight in [1.0, 2.0]
48
      @testset "primal_weight=$primal_weight" begin
49
        params = PrimalDualOptimizerParameters(
50
          RecoverPrimalDualMethodFromString(method_name), # method
51
          0.2,  # step_size
52
          primal_weight,  # primal_weight
53
          1, # record every
54
          20,   # printing frequency
55
          false, # verbose
56
          2000,  # iteration limit
57
          NoRestarts(), # restart scheme
58
          nothing,
59
          1e-8,
60
        )
61

62
        ##############
63
        # simple lp1 #
64
        ##############
65
        # No restarts
66
        solver_output = optimize(params, simple_lp1, [1.0], [1.0])
67
        @test norm(solver_output.primal_solution) < 1e-7
68
        @test norm(solver_output.dual_solution) < 1e-7
69
        @test solver_output.status == STATUS_OPTIMAL
70
        @test solver_output.iteration_stats.iteration[end] < 1500
71

72
        # Fixed frequency restarts
73
        params.restart_scheme = FixedFrequencyRestarts(restart_length, true)
74
        solver_output = optimize(params, simple_lp1, [1.0], [1.0])
75
        @test norm(solver_output.primal_solution) < 1e-7
76
        @test norm(solver_output.dual_solution) < 1e-7
77
        @test solver_output.status == STATUS_OPTIMAL
78
        @test solver_output.iteration_stats.iteration[end] < 600
79

80
        # Adaptive restarts
81
        params.restart_scheme = AdaptiveRestarts(exp(-1), true)
82
        solver_output = optimize(params, simple_lp1, [1.0], [1.0])
83
        @test norm(solver_output.primal_solution) < 1e-7
84
        @test norm(solver_output.dual_solution) < 1e-7
85
        @test solver_output.status == STATUS_OPTIMAL
86
        @test solver_output.iteration_stats.iteration[end] < 600
87

88
        # Adaptive restarts
89
        params.restart_scheme = AdaptiveRestarts(exp(-1), false)
90
        solver_output = optimize(params, simple_lp1, [1.0], [1.0])
91
        @test norm(solver_output.primal_solution) < 1e-7
92
        @test norm(solver_output.dual_solution) < 1e-7
93
        @test solver_output.status == STATUS_OPTIMAL
94
        @test solver_output.iteration_stats.iteration[end] < 600
95

96
        ##############
97
        # simple lp2 #
98
        ##############
99
        # No restarts
100
        solver_output = optimize(params, simple_lp2)
101
        @test norm(solver_output.primal_solution - [0.0, 1.0]) < 1e-3
102
        @test norm(solver_output.dual_solution - [1.0]) < 1e-3
103
        @test solver_output.status == STATUS_OPTIMAL
104

105
        # Fixed frequency restarts
106
        params.restart_scheme = FixedFrequencyRestarts(restart_length, true)
107
        solver_output = optimize(params, simple_lp2)
108
        @test norm(solver_output.primal_solution - [0.0, 1.0]) < 1e-7
109
        @test norm(solver_output.dual_solution - [1.0]) < 1e-7
110
        @test solver_output.status == STATUS_OPTIMAL
111

112
        # Adaptive restarts
113
        params.restart_scheme = AdaptiveRestarts(exp(-1), true)
114
        solver_output = optimize(params, simple_lp2)
115
        @test norm(solver_output.primal_solution - [0.0, 1.0]) < 1e-7
116
        @test norm(solver_output.dual_solution - [1.0]) < 1e-7
117
        @test solver_output.status == STATUS_OPTIMAL
118

119
        #######
120
        # lp1 #
121
        #######
122
        # No restarts
123
        params.kkt_tolerance = 1e-6
124
        params.restart_scheme = NoRestarts()
125
        solver_output = optimize(params, lp1)
126
        @test norm(
127
          solver_output.primal_solution - lp1_primal_optimal_solution,
128
        ) < 1e-2
129
        @test norm(solver_output.dual_solution - lp1_dual_optimal_solution) <
130
              1e-2
131
        # PDHG with NoRestarts hits the termination limit, so we have a
132
        # larger tolerance on the primal and dual solutions, and don't test
133
        # status.
134

135
        # Fixed frequency restarts
136
        params.restart_scheme = FixedFrequencyRestarts(restart_length, true)
137
        solver_output = optimize(params, lp1)
138
        @test norm(
139
          solver_output.primal_solution - lp1_primal_optimal_solution,
140
        ) < 1e-5
141
        @test norm(solver_output.dual_solution - lp1_dual_optimal_solution) <
142
              1e-5
143
        @test solver_output.status == STATUS_OPTIMAL
144

145
        # Adaptive restarts
146
        params.restart_scheme = AdaptiveRestarts(exp(-1), true)
147
        solver_output = optimize(params, lp1)
148
        @test norm(
149
          solver_output.primal_solution - lp1_primal_optimal_solution,
150
        ) < 1e-5
151
        @test norm(solver_output.dual_solution - lp1_dual_optimal_solution) <
152
              1e-5
153
        @test solver_output.status == STATUS_OPTIMAL
154
      end
155
    end
156
  end
157
end
158