tutorial: warm-start + sparse inverse covariance + MOSEK appendix; split bootstrap; renumber

Content:
- ch8 (DPP): add "warm starts" section — warm-start cuts solver time,
  DPP cuts compile time; orthogonal and complementary, with a runnable
  warm-start sweep verified to match the cold path.
- ch9 (portfolio): make the efficient-frontier sweep a real Parameter
  sweep (build once, set value(gamma)) instead of rebuilding each gamma.
- new ch10: Sparse Inverse Covariance Estimation (log_det + l1 budget,
  PSD variable), verified vs CVXR 1.9.1; base chol() sampling, no expm.
- split the old combined "Showcase" chapter into ch12 Image
  Reconstruction (TV in-painting only) and ch13 The Bootstrap
  (expanded — general inference recipe for any psolve() estimator).
- new Appendix A: Using a Commercial Solver (MOSEK) — get/use it plus a
  static, pre-rendered verbose example (genuine MOSEK 11.1.2 SDP log).

Structure:
- renumber all chapter files to sequential integers (drop 09a/10a/11a
  suffixes); appendix-mosek.qmd; References now follows the appendix.
- code-overflow: scroll so long code lines get a horizontal scrollbar.

Author: bnaras

08-dpp.qmd

@@ -154,6 +154,44 @@ dim(betas)   # m coefficients x number of lambda values
 Whenever you plan to solve the *same shape* of problem many times --- a path, a
 grid search, a simulation, a bootstrap --- reach for a `Parameter`.
 
+## The other half of the solve: warm starts
+
+Recall that every solve splits into a **compile time** and a **solver time**. DPP
+attacks the *compile* time --- it reuses the canonicalized problem structure so we
+never re-stuff the matrices. But the solver still starts **cold** each time, refining
+an internal guess from scratch until it converges.
+
+A **warm start** attacks the *other* half. A numerical solver is iterative; if you
+hand it a starting point already near the answer, it converges in fewer iterations.
+In a sweep over a fine grid the solution at one $\lambda$ is an excellent guess for
+the next, so the previous answer is exactly the start point we want. In CVXR that is
+one flag --- `warm_start = TRUE` uses the current variable values as the solver's
+initial point:
+
+```{r}
+#| label: warm-sweep
+betas_warm <- sapply(lambda_vals, function(l) {
+  value(lambda) <- l
+  psolve(path_prob, warm_start = TRUE)   # start from the previous lambda's solution
+  value(beta)
+})
+max(abs(betas_warm - betas))             # same path, to solver tolerance
+```
+
+The two ideas are **orthogonal and complementary**, not the same trick:
+
+| | reuses | cuts | wins biggest when |
+|---|---|---|---|
+| **DPP** | the compiled problem structure | *compile* time | compilation dominates --- small problems solved *many* times |
+| **Warm start** | the previous solution | *solver* time | the solver dominates --- large/expensive problems, close-together solves |
+
+DPP does nothing for the solver's iterations; a warm start does nothing for
+compilation. Each pays off precisely where the other does not, so a tight sweep wants
+**both** --- compile once (DPP), then step the solver from each previous answer
+(warm start). Warm-starting is solver-dependent (about 9 of CVXR's 15 solvers support
+it --- OSQP, SCS, Clarabel, MOSEK, Gurobi, PIQP, HiGHS, and more); a solver that
+cannot warm-start simply ignores the flag.
+
 ## Exercise
 
 **Exercise (⭐).** The expression `gamma * p_norm(x, 1)` is DPP, but
@@ -182,5 +220,7 @@ parameter should be computed in R and stored in another parameter.
   CVXR compiles once and re-solves cheaply.
 - Use it for any **repeated solve** --- paths, grid searches, simulations. It is the
   efficient way to do the sweep you met in @sec-elastic-net.
+- **DPP cuts compile time; a warm start cuts solver time.** They are independent ---
+  combine `Parameter`s with `warm_start = TRUE` to speed up both halves of a sweep.
 
 For more, see the [DPP Tutorial](https://cvxr.rbind.io/examples/dpp/tutorial.html).

09-portfolio.qmd

@@ -26,6 +26,7 @@ underlying models are due to @Markowitz:1952 (mean-variance) and
 #| include: false
 library(CVXR)
 library(ggplot2)
+library(RColorBrewer)
 source("_check_status.R")
 ```
 
@@ -63,19 +64,23 @@ colnames(R)
 ## Trace the efficient frontier
 
 For each $\gamma$ we solve the Markowitz problem and record the realized risk and
-return. As in @sec-elastic-net, this is the **same problem re-solved** over a grid ---
-exactly the place where the `Parameter` trick from @sec-dpp pays off.
+return. This is the **same problem re-solved** over a grid --- exactly where the
+`Parameter` trick from @sec-dpp pays off. We make $\gamma$ a `Parameter`, build the
+problem **once**, and each iteration just sets `value(gamma)` and re-solves, reusing
+the cached compilation instead of canonicalizing from scratch 40 times.
 
 ```{r}
 #| label: frontier
-w    <- Variable(n)
-ret  <- t(mu) %*% w
-risk <- quad_form(w, Sigma)               # w' Sigma w, a convex quadratic
-cons <- list(w >= 0, sum(w) == 1)
+w     <- Variable(n)
+ret   <- t(mu) %*% w
+risk  <- quad_form(w, Sigma)              # w' Sigma w, a convex quadratic
+cons  <- list(w >= 0, sum(w) == 1)
+gamma <- Parameter(nonneg = TRUE)         # the risk-aversion knob, set per solve
+prob  <- Problem(Maximize(ret - gamma * risk), cons)   # built once
 
 gammas <- 10^seq(-1, 3, length.out = 40)
 fr <- sapply(gammas, function(g) {
-  prob <- Problem(Maximize(ret - g * risk), cons)
+  value(gamma) <- g                       # swap in the new gamma; no rebuild
   psolve(prob)
   c(risk = sqrt(value(risk)), ret = value(ret), w = value(w))
 })
@@ -101,6 +106,36 @@ ggplot(df, aes(risk, ret)) +
 Every point on the blue curve is the best return achievable at that level of risk; the
 individual stocks (red) all sit *below* it --- diversification at work.
 
+The frontier shows the risk-return *outcome*, but not *what the portfolio holds*. Since
+`fr` already captured the weight vector at every $\gamma$, we can read the composition
+straight back out and watch diversification happen. At a few points along the grid we
+stack each asset's share of the budget:
+
+```{r}
+#| label: frontier-composition
+#| fig-cap: "Portfolio composition along the frontier. Small gamma concentrates the budget in a few assets; large gamma spreads it out --- diversification, made visible."
+markers <- c(1, 10, 20, 30, 40)
+wmat <- fr[grep("^w", rownames(fr)), markers]   # 12 assets x 5 values of gamma
+rownames(wmat) <- names(mu)
+
+comp <- data.frame(
+  gamma  = factor(rep(signif(gammas[markers], 2), each = n),
+                  levels = signif(gammas[markers], 2)),
+  asset  = factor(rep(names(mu), times = length(markers)), levels = names(mu)),
+  weight = as.vector(wmat))
+
+ggplot(comp, aes(gamma, weight, fill = asset)) +
+  geom_col() +
+  scale_fill_manual(values = brewer.pal(n, "Paired")) +
+  labs(x = expression("Risk aversion  " * gamma),
+       y = "Fraction of budget", fill = "Asset") +
+  theme_minimal()
+```
+
+The left-most bar (greedy, low $\gamma$) leans on one or two names; as $\gamma$ grows the
+bars fragment into the full palette. That shift *is* the efficient frontier, seen from
+the portfolio's side.
+
 ## A different risk measure: CVaR
 
 Variance penalizes upside and downside alike. **Conditional Value-at-Risk** (CVaR)
@@ -140,6 +175,38 @@ Same data, same `CVXR` machinery, a different notion of risk --- and a different
 typically more tail-aware, allocation. Switching risk measures is, once again, just a
 change to the objective.
 
+To see *how* the allocation differs, hold the return target fixed and put the two
+side by side. The fair mean-variance counterpart minimizes variance subject to the same
+return floor, so both bars sit at the same expected return and only the risk measure
+changes:
+
+```{r}
+#| label: cvar-vs-mv
+#| fig-cap: "Same return target, two risk measures. Variance and CVaR reallocate the budget differently."
+## mean-variance allocation at the SAME return target, for an apples-to-apples comparison
+wv_mv   <- Variable(n)
+prob_mv <- Problem(Minimize(quad_form(wv_mv, Sigma)),
+                   list(sum(wv_mv) == 1, wv_mv >= 0, t(mu) %*% wv_mv >= target))
+psolve(prob_mv)
+check_solver_status(prob_mv)
+
+comp2 <- data.frame(
+  measure = factor(rep(c("Mean-variance", "Mean-CVaR"), each = n),
+                   levels = c("Mean-variance", "Mean-CVaR")),
+  asset   = factor(rep(names(mu), 2), levels = names(mu)),
+  weight  = c(value(wv_mv), value(wv)))
+
+ggplot(comp2, aes(measure, weight, fill = asset)) +
+  geom_col() +
+  scale_fill_manual(values = brewer.pal(n, "Paired")) +
+  labs(x = NULL, y = "Fraction of budget", fill = "Asset") +
+  theme_minimal()
+```
+
+Both portfolios earn the same expected return, yet CVaR --- caring only about the bad
+tail --- tilts the budget toward a different mix than variance, which penalizes upside
+and downside alike.
+
 ## Exercise
 
 **Exercise (⭐).** Add a **diversification** rule to the mean-variance problem: no

10-invcov.qmd

@@ -0,0 +1,141 @@
+# Sparse Inverse Covariance Estimation {#sec-invcov}
+
+<!-- LIVE. Adapted from cvxr_docs statistics/sparse-inverse-covariance-estimation.qmd
+     (Fu & Narasimhan). Sampling uses base chol() instead of expm::sqrtm to keep
+     dependencies light. Verified vs CVXR 1.9.1 (2026-06-10): true nonzeros = 24;
+     alpha=10 -> 46, alpha=4 -> 24 (recovers truth), alpha=1 -> 12. -->
+
+::: callout-note
+## What you'll be able to do
+Estimate the **graph** of conditional dependencies among variables --- which ones are
+related once you account for all the others --- by fitting a **sparse inverse
+covariance** matrix with a log-determinant objective and an $\ell_1$ budget.
+:::
+
+```{r}
+#| label: setup
+#| include: false
+library(CVXR)
+library(Matrix)
+library(ggplot2)
+source("_check_status.R")
+```
+
+## The idea
+
+The portfolio chapter (@sec-portfolio) took a covariance matrix as given. Here we go
+the other way and **estimate** it --- and specifically its inverse, which carries the
+more interesting story.
+
+For Gaussian data, the **inverse** covariance $S = \Sigma^{-1}$ (the *precision*
+matrix) encodes **conditional independence**: $S_{ij} = 0$ exactly when variables $i$
+and $j$ are independent *given all the others*. So a **sparse** $S$ is a sparse graph
+of direct relationships --- often what we actually want to discover. With more
+variables than the data can pin down, we impose that sparsity directly.
+
+Given i.i.d. samples $x_i \sim N(0, \Sigma)$ with sample covariance $Q$, we maximize
+the Gaussian log-likelihood while capping the total size of $S$:
+
+$$
+\begin{array}{ll}
+\underset{S}{\mbox{maximize}} & \log\det(S) - \mathop{\mathrm{tr}}(SQ) \\
+\mbox{subject to} & S \succeq 0, \qquad \sum_{i,j} |S_{ij}| \le \alpha .
+\end{array}
+$$
+
+The $\ell_1$ budget $\alpha \ge 0$ controls sparsity: smaller $\alpha$ forces more
+off-diagonal entries to exactly zero. The objective is concave (`log_det` is concave;
+the trace term is linear) and the constraints are convex, so CVXR solves it directly.
+
+## Synthetic data with a known answer
+
+So we can check the estimate, we build a *true* sparse precision matrix $S_{\text{true}}$,
+invert it to get the covariance $R$, and sample from $N(0, R)$:
+
+```{r}
+#| label: data
+set.seed(1)
+n <- 10        # number of variables
+m <- 1000      # number of samples
+
+A      <- rsparsematrix(n, n, 0.15, rand.x = rnorm)
+Strue  <- as.matrix(A %*% t(A) + 0.05 * diag(n))   # sparse, symmetric, positive definite
+R      <- solve(Strue)                              # the covariance to sample from
+
+x <- matrix(rnorm(m * n), m, n) %*% chol(R)         # rows are i.i.d. N(0, R)
+Q <- cov(x)                                          # sample covariance
+
+sum(abs(Strue) > 1e-4)                               # how many true nonzeros (of 100)
+```
+
+## Fit at several sparsity budgets
+
+A semidefinite variable is declared with `PSD = TRUE`. We sweep $\alpha$ and, after
+each solve, zero out numerically negligible entries to read off the recovered pattern:
+
+```{r}
+#| label: fit
+alphas <- c(10, 6, 4, 1)
+S      <- Variable(c(n, n), PSD = TRUE)
+obj    <- Maximize(log_det(S) - matrix_trace(S %*% Q))
+
+S_est <- lapply(alphas, function(a) {
+  prob <- Problem(obj, list(sum(abs(S)) <= a))
+  psolve(prob)
+  check_solver_status(prob)
+  Sa <- value(S)
+  Sa[abs(Sa) <= 1e-4] <- 0            # treat tiny entries as exact zeros
+  Sa
+})
+sapply(S_est, function(M) sum(M != 0))  # nonzeros recovered at each alpha
+```
+
+Note `matrix_trace(S %*% Q)` rather than `sum(diag(S %*% Q))` --- same value, but the
+dedicated atom keeps the expression tree compact. And `log_det(S)` is its own atom:
+there is no `det()` for a CVXR variable, because the log-determinant is the convex
+object DCP understands.
+
+## What the budget buys
+
+The sparsity pattern of the estimate tightens as $\alpha$ shrinks. We plot the nonzero
+mask of the truth alongside each fit:
+
+```{r}
+#| label: pattern-plot
+#| fig-width: 9
+#| fig-height: 2.4
+#| fig-cap: "Nonzero pattern (dark = nonzero) of the true precision matrix and the estimates. Smaller alpha -> sparser. Around alpha = 4 the estimate matches the truth's sparsity."
+mask_df <- function(M, lab) {
+  data.frame(i = rep(seq_len(n), n), j = rep(seq_len(n), each = n),
+             nz = factor(as.integer(as.vector(M) != 0)), panel = lab)
+}
+panels <- rbind(
+  mask_df(Strue, "truth"),
+  do.call(rbind, Map(function(M, a) mask_df(M, paste0("alpha = ", a)), S_est, alphas)))
+panels$panel <- factor(panels$panel, levels = c("truth", paste0("alpha = ", alphas)))
+
+ggplot(panels, aes(i, j, fill = nz)) +
+  geom_tile(color = "grey80") +
+  scale_fill_manual(values = c("0" = "white", "1" = "steelblue"), guide = "none") +
+  scale_y_reverse() + coord_equal() +
+  facet_wrap(~ panel, nrow = 1) +
+  theme_void() + theme(strip.text = element_text(margin = margin(4, 0, 4, 0)))
+```
+
+At $\alpha = 10$ the budget is loose and the estimate is dense; as $\alpha$ falls the
+off-diagonal entries are squeezed to zero, and near $\alpha = 4$ the recovered pattern
+matches the true number of nonzeros. Too small ($\alpha = 1$) and we lose real edges.
+In practice you would pick $\alpha$ by cross-validation --- and because the only thing
+changing across the sweep is one constant, this is a textbook case for the `Parameter`
+trick of @sec-dpp.
+
+## Takeaways
+
+- The **inverse** covariance is the object of interest: its zeros are **conditional
+  independencies** --- a graph of direct relationships.
+- `log_det(S) - matrix_trace(S %*% Q)` with `S` declared `PSD = TRUE` and an
+  $\ell_1$ budget is the whole model --- a few lines of CVXR.
+- One knob, $\alpha$, trades density for sparsity; sweep it with a `Parameter`.
+
+See [Sparse Inverse Covariance Estimation](https://cvxr.rbind.io/examples/statistics/sparse-inverse-covariance-estimation.html)
+for the full treatment.