Вот как это можно сделать с CVXR:
> library(CVXR)
>
> #
> # data
> #
> Y <- c(167, 136, 195, 174, 144, 135, 89, 81, 114, 111)
> X1 <- c(2.9, 3.4, 0.7, 1.1, 3.5, 5.0, 6.7, 4.7, 3.7, 8.8)
> X2 <- c(60, 47, 63, 62, 40, 60, 50, 35, 40, 40)
>
> #
> # organize as matrix
> #
> X <- cbind(1,X1,X2)
> X
X1 X2
[1,] 1 2.9 60
[2,] 1 3.4 47
[3,] 1 0.7 63
[4,] 1 1.1 62
[5,] 1 3.5 40
[6,] 1 5.0 60
[7,] 1 6.7 50
[8,] 1 4.7 35
[9,] 1 3.7 40
[10,] 1 8.8 40
>
>
> #
> # standard regression
> #
> lm(Y~X1+X2)
Call:
lm(formula = Y ~ X1 + X2)
Coefficients:
(Intercept) X1 X2
82.616 -7.716 1.675
>
> #
> # standard regression as QP
> #
> beta <- Variable(3)
>
> problem <- Problem(Minimize(sum_squares(Y - X %*% beta)))
> result <- solve(problem, verbose=T)
-----------------------------------------------------------------
OSQP v0.6.0 - Operator Splitting QP Solver
(c) Bartolomeo Stellato, Goran Banjac
University of Oxford - Stanford University 2019
-----------------------------------------------------------------
problem: variables n = 13, constraints m = 10
nnz(P) + nnz(A) = 50
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 1.00e-01 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 0.0000e+00 1.95e+02 1.70e+05 1.00e-01 3.47e-04s
50 2.8959e+03 3.20e-07 2.29e-06 1.00e-01 8.43e-04s
plsh 2.8959e+03 4.21e-14 4.51e-13 --------- 1.36e-03s
status: solved
solution polish: successful
number of iterations: 50
optimal objective: 2895.8618
run time: 1.36e-03s
optimal rho estimate: 4.57e-02
> result$getValue(beta)
[,1]
[1,] 82.616486
[2,] -7.716095
[3,] 1.674722
>
> #
> # add constraints
> # beta[1] >= 0 (constant term)
> # beta[2] >= 0.0001 (instead of > we do >= with a small tolerance)
> # beta[3] <= 0.9999 (instead of < we do <= with a small tolerance)
> #
> problem <- Problem(Minimize(sum_squares(Y - X %*% beta)),
+ list(beta[1]>=0,
+ beta[2]>=0.0001,
+ beta[3]<=0.9999))
> result <- solve(problem, verbose=T)
-----------------------------------------------------------------
OSQP v0.6.0 - Operator Splitting QP Solver
(c) Bartolomeo Stellato, Goran Banjac
University of Oxford - Stanford University 2019
-----------------------------------------------------------------
problem: variables n = 13, constraints m = 13
nnz(P) + nnz(A) = 53
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 1.00e-01 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 0.0000e+00 1.95e+02 1.71e+05 1.00e-01 5.30e-04s
150 7.8534e+03 9.37e-05 5.05e-04 4.62e+00 1.51e-03s
plsh 7.8533e+03 5.55e-14 2.84e-14 --------- 1.94e-03s
status: solved
solution polish: successful
number of iterations: 150
optimal objective: 7853.3365
run time: 1.94e-03s
optimal rho estimate: 2.98e+00
> result$getValue(beta)
[,1]
[1,] 84.90456
[2,] 0.00010
[3,] 0.99990
>