2222# ```math
2323# \begin{split}
2424# \begin{array} {ll}
25- # \mbox{minimize} & e^{\top}e + \alpha (w^2 + b^2 ) \\
25+ # \mbox{minimize} & e^{\top}e + \alpha (w^2) \\
2626# \mbox{s.t.} & e_{i} = y_{i} - w x_{i} - b \quad \quad i=1..N \\
2727# \end{array}
2828# \end{split}
@@ -36,15 +36,15 @@ import DiffOpt
3636import Random
3737import Ipopt
3838import Plots
39- import LinearAlgebra: normalize!, dot
39+ using LinearAlgebra: dot
4040
4141# ## Define and solve the problem
4242
4343# Construct a set of noisy (guassian) data points around a line.
4444
4545Random. seed! (42 )
4646
47- N = 100
47+ N = 150
4848
4949w = 2 * abs (randn ())
5050b = rand ()
@@ -64,76 +64,106 @@ function fit_ridge(X, Y, alpha = 0.1)
6464 set_silent (model)
6565 @variable (model, w) # angular coefficient
6666 @variable (model, b) # linear coefficient
67- @variable (model, e[1 : N]) # approximation error
68- # # constraint defining approximation error
69- @constraint (model, cons[i= 1 : N], e[i] == Y[i] - w * X[i] - b)
67+ # # expression defining approximation error
68+ @expression (model, e[i= 1 : N], Y[i] - w * X[i] - b)
7069 # # objective minimizing squared error and ridge penalty
7170 @objective (
7271 model,
7372 Min,
74- dot (e, e) + alpha * (sum (w * w) + sum (b * b) ),
73+ 1 / N * dot (e, e) + alpha * (w ^ 2 ),
7574 )
7675 optimize! (model)
77- return model, w, b, cons # return model, variables and constraints references
76+ return model, w, b # return model & variables
7877end
7978
8079
81- # Train on the data generated.
80+ # Plot the data points and the fitted line for different alpha values
8281
83- model, w, b, cons = fit_ridge (X, Y)
84- ŵ, b̂ = value (w), value (b)
85-
86- # We can visualize the approximating line.
87-
88- p = Plots. scatter (X, Y, label= " "
82+ p = Plots. scatter (X, Y, label= nothing , legend= :topleft )
8983mi, ma = minimum (X), maximum (X)
90- Plots. plot! (p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :red , label = " "
84+ Plots. title! ( " Fitted lines and points "
9185
86+ for alpha in 0.5 : 0.5 : 1.5
87+ local model, w, b = fit_ridge (X, Y, alpha)
88+ ŵ = value (w)
89+ b̂ = value (b)
90+ Plots. plot! (p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], label= " alpha=$alpha " = 2 )
91+ end
92+ p
9293
9394# ## Differentiate
9495
9596# Now that we've solved the problem, we can compute the sensitivity of optimal
96- # values of the angular coefficient `w` with
97+ # values of the slope `w` with
9798# respect to perturbations in the data points (`x`,`y`).
9899
99- # Begin differentiating the model.
100- # analogous to varying θ in the expression:
101- # ```math
102- # e_i = (y_{i} + \theta_{y_i}) - w (x_{i} + \theta_{x_{i}}) - b
103- # ```
100+ alpha = 0.4
101+ model, w, b = fit_ridge (X, Y, alpha)
102+ ŵ = value (w)
103+ b̂ = value (b)
104+
105+ # We first compute sensitivity of the slope with respect to a perturbation of the independent
106+ # variable `x`.
104107
105- ∇ = zero (X)
108+ # Recalling that the points $(x_i, y_i)$ appear in the objective function as:
109+ # `(yi - b - w*xi)^2`, the `DiffOpt.ForwardInObjective` attribute must be set accordingly,
110+ # with the terms multiplying the parameter in the objective.
111+ # When considering the perturbation of a parameter θ, `DiffOpt.ForwardInObjective()` takes in the expression in the
112+ # objective that multiplies θ.
113+ # If θ appears with a quadratic and a linear form: `θ^2 a x + θ b y`, then the expression to pass to
114+ # `ForwardInObjective` is `2θ a x + b y`.
115+
116+ # Sensitivity with respect to x and y
117+
118+ ∇y = zero (X)
119+ ∇x = zero (X)
106120for i in 1 : N
107- for j in 1 : N
108- MOI. set (
109- model,
110- DiffOpt. ForwardInConstraint (),
111- cons[j],
112- i == j ? index (w) + 1.0 : 0.0 * index (w)
113- )
114- end
121+ MOI. set (
122+ model,
123+ DiffOpt. ForwardInObjective (),
124+ 2 w^ 2 * X[i] + 2 b * w - 2 * w * Y[i]
125+ )
126+ DiffOpt. forward (model)
127+ ∇x[i] = MOI. get (
128+ model,
129+ DiffOpt. ForwardOutVariablePrimal (),
130+ w
131+ )
132+ MOI. set (
133+ model,
134+ DiffOpt. ForwardInObjective (),
135+ (2 Y[i] - 2 b - 2 w * X[i]),
136+ )
115137 DiffOpt. forward (model)
116- dw = MOI. get (
138+ ∇y[i] = MOI. get (
117139 model,
118140 DiffOpt. ForwardOutVariablePrimal (),
119141 w
120142 )
121- ∇[i] = abs (dw)
122143end
123144
124- normalize! (∇);
145+ # Visualize point sensitivities with respect to regression points.
146+
147+ p = Plots. scatter (
148+ X, Y,
149+ color = [dw < 0 ? :blue : :red for dw in ∇x],
150+ markersize = [5 * abs (dw) + 1.2 for dw in ∇x],
151+ label = " "
152+ )
153+ mi, ma = minimum (X), maximum (X)
154+ Plots. plot! (p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :blue , label = " "
155+ Plots. title! (" Regression slope sensitivity with respect to x"
125156
126- # Visualize point sensitivities with respect to regressing line.
127- # Note that the gradients are normalized.
157+ #
128158
129159p = Plots. scatter (
130160 X, Y,
131- color = [x > 0 ? :red : :blue for x in ∇],
132- markersize = [25 * abs (x) for x in ∇],
161+ color = [dw < 0 ? :blue : :red for dw in ∇y ],
162+ markersize = [5 * abs (dw) + 1.2 for dw in ∇y ],
133163 label = " "
134164)
135165mi, ma = minimum (X), maximum (X)
136- Plots. plot! (p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :red , label = " "
166+ Plots. plot! (p, [mi, ma], [mi * ŵ + b̂, ma * ŵ + b̂], color = :blue , label = " "
167+ Plots. title! (" Regression slope sensitivity with respect to y"
137168
138- # Note the points in the extremes of the line segment are larger because
139- # moving those points has a stronger effect on the angular coefficient of the line.
169+ # Note the points with less central `x` values induce a greater y sensitivity of the slope.
0 commit comments