Mean-variance optimization

Let there be n assets. Collect their expected returns in a vector μ (length n), their volatilities in a vector σ, and their pairwise correlations in a matrix ρ (ones on the diagonal, symmetric). The covariance matrix is built exactly the way M3.02 §4.3 built a single covariance, just for every pair at once:

Σ = diag(σ) · ρ · diag(σ)          (so that Σ_ij = σ_i σ_j ρ_ij, and Σ_ii = σ_i²)

A portfolio is a weight vector w (length n, summing to 1). Its expected return and variance are the direct matrix generalization of formulas you already own:

Portfolio return:    Rp = w′μ = Σ_i w_i μ_i                              (a weighted average — same as M3.02)
Portfolio variance:  σp² = w′Σw = Σ_i Σ_j w_i w_j σ_i σ_j ρ_ij           (the double sum M3.02 §4.3–4.4 wrote out by hand)

How is mean-variance optimization solved?

Expand w′Σw for two assets and you get back M3.02's exact two-asset formula, this is not a new idea, it is the old one written so it scales. The full Markowitz problem is: minimize w′Σw subject to w′μ = target return and w′1 = 1 (and, when the mandate demands it, bounds on individual weights: see constrained mean-variance optimization). This is a quadratic program. It has a closed form when the only constraints are the two equalities above (solved via a bordered linear system, a Lagrangian you can write out but need not memorize), and it needs a numerical solver: scipy.optimize.minimize, the instant you add a bound like "no shorting."

Modern Data Analysis, the Python way. Every frontier, tangency portfolio, and sensitivity number in this node came from a script no longer than this: build Sigma = np.diag(sigma) @ corr @ np.diag(sigma); get the global-minimum-variance weights from Sinv = np.linalg.inv(Sigma); w = Sinv @ ones / (ones @ Sinv @ ones); get an unconstrained frontier point at a target return by solving the bordered KKT system with np.linalg.solve; get the tangency portfolio from w = Sinv @ (mu - rf) renormalized to sum to 1; and get any constrained version, long-only, capped, whatever the mandate requires, by handing scipy.optimize.minimize(lambda w: w @ Sigma @ w, ..., method="SLSQP", bounds=..., constraints=...) the same objective with bounds bolted on. That is the entire mechanical core of institutional portfolio-construction software, a dozen lines, not a black box. Two things worth checking every single time, because no error message will catch either for you: confirm your correlation matrix is a valid input at all (np.linalg.eigvalsh(corr), every eigenvalue must be positive, see assembling the machine), and confirm a GMV solution's first-order condition (Sigma @ w should be a constant vector) before trusting anything built on top of it. The full production version of this exact code, real constraint sets, robust and shrinkage estimators, Black-Litterman, risk parity, and the resampling method that directly answers the fragility shown in the input-sensitivity failure: is QD1.04 (quant-code region); this callout shows you the shape of it so that module reads as an extension of something you already understand, not a new subject.