1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
---
title: "ECOSolveR Examples"
date: '`r Sys.Date()`'
output:
html_document:
fig_caption: yes
theme: cerulean
toc: yes
toc_depth: 2
vignette: >
%\VignetteIndexEntry{ECOSolveR Examples}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r echo=F}
### get knitr just the way we like it
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE,
error = FALSE,
tidy = FALSE,
cache = FALSE
)
```
## $L_1$ minimization (Linear Programming)
We solve the following problem that arises for example in sparse
signal reconstruction problems such as compressed sensing:
$$
\mbox{minimize } ||x||_1 \mbox{ ($L_1$) }\\
\mbox{subject to } Ax = b
$$
with $x\in R^n$, $A \in R^{m \times n}$ and $m\leq n.$ Reformulate the
problem expressing the $L_1$ norm of $x$ as follows
$$
x \leq u\\
-x \leq u\\
$$
where $u\in R^n$ and we minimize the sum of $u$. The reformulated
problem using the stacked variables
$$
z = \begin{pmatrix}x\\u\end{pmatrix}
$$
is now
$$
\mbox{minimize } c^{\top}z\\
\mbox{subject to } \tilde{A}x = b \mbox{ (LP) }\\
Gx \leq h
$$
where the inequality is with respective to the positive orthant.
Here is the R code that generates a random instance of this problem
and solves it.
```{r}
library(ECOSolveR)
library(Matrix)
set.seed(182391)
n <- 1000L
m <- 10L
density <- 0.01
c <- c(rep(0.0, n), rep(1.0, n))
```
First, a function to generate random sparse matrices with normal
entries.
```{r}
sprandn <- function(nrow, ncol, density) {
items <- ceiling(nrow * ncol * density)
matrix(c(rnorm(items),
rep(0, nrow * ncol - items)),
nrow = nrow)
}
```
```{r}
A <- sprandn(m, n, density)
Atilde <- Matrix(cbind(A, matrix(rep(0.0, m * n), nrow = m)), sparse = TRUE)
b <- rnorm(m)
I <- diag(n)
G <- rbind(cbind(I, -I),
cbind(-I, -I))
G <- as(G, "dgCMatrix")
h <- rep(0.0, 2L * n)
dims <- list(l = 2L * n, q = NULL, e = 0L)
```
Note how ECOS expects sparse matrices, not ordinary matrices.
```{r}
## Solve the problem
z <- ECOS_csolve(c = c, G = G, h = h, dims = dims, A = Atilde, b = b)
```
We check that the solution was found.
```{r}
names(z)
z$infostring
```
Extract the solution.
```{r}
x <- z$x[1:n]
u <- z$x[(n+1):(2*n)]
nnzx = sum(abs(x) > 1e-8)
sprintf("x reconstructed with %d non-zero entries", nnzx / length(x) * 100)
```
## Parametric Optimization with the Update-Solve Loop
When solving a sequence of related problems that share the same
sparsity structure but differ in numerical data, the multi-step
lifecycle API avoids repeating the expensive symbolic analysis and
KKT ordering phase. The pattern is:
1. `ECOS_setup()` — one-time setup (symbolic analysis)
2. Loop: `ECOS_update()` + `ECOS_solve()` — cheap re-solves
3. `ECOS_cleanup()` — free workspace
Here we vary the right-hand side $h$ of the LP from the previous
example, tracing how the optimal value changes as we scale $h$.
```{r}
## Set up workspace once
ws <- ECOS_setup(c = c, G = G, h = h, dims = dims, A = Atilde, b = b)
## Sweep a parameter: scale h from 0 to 1
alphas <- seq(0, 1, length.out = 11)
pcosts <- numeric(length(alphas))
for (i in seq_along(alphas)) {
ECOS_update(ws, h = h + alphas[i])
res <- ECOS_solve(ws)
pcosts[i] <- res$summary[["pcost"]]
}
ECOS_cleanup(ws)
## Show results
data.frame(alpha = alphas, pcost = round(pcosts, 4))
```
For comparison, the same sweep using `ECOS_csolve()` in a loop
would call `ECOS_setup` and `ECOS_cleanup` at every iteration.
The lifecycle API is particularly beneficial when the problem
dimensions are large and setup dominates solve time.
|
