From e8191a64dd430e3351b99c9906d4e828666121f6 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 10:56:08 +1000 Subject: [PATCH 1/9] doc link --- lectures/linear_algebra.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 33f3dc53e..2cfa0d597 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -65,7 +65,7 @@ These are the kinds of topics addressed by linear algebra. In this lecture we will cover the basics of linear and matrix algebra, treating both theory and computation. -We admit some overlap with [this lecture](https://python-programming.quantecon.org/numpy.html), where operations on NumPy arrays were first explained. +We admit some overlap with {doc}`programming:numpy`, where operations on NumPy arrays were first explained. Note that this lecture is more theoretical than most, and contains background material that will be used in applications as we go along. From ff470355006899d592ce039c9121c437ace79f6a Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 11:15:58 +1000 Subject: [PATCH 2/9] Revert "doc link" This reverts commit c45a974d503beb7e502223275d607d4565cbe846. --- lectures/linear_algebra.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 2cfa0d597..33f3dc53e 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -65,7 +65,7 @@ These are the kinds of topics addressed by linear algebra. In this lecture we will cover the basics of linear and matrix algebra, treating both theory and computation. -We admit some overlap with {doc}`programming:numpy`, where operations on NumPy arrays were first explained. +We admit some overlap with [this lecture](https://python-programming.quantecon.org/numpy.html), where operations on NumPy arrays were first explained. Note that this lecture is more theoretical than most, and contains background material that will be used in applications as we go along. From 9821db79b6ea456032c3fed9a4987209164c8e69 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 11:37:57 +1000 Subject: [PATCH 3/9] update figure style --- lectures/linear_algebra.md | 41 ++++++++++++++++---------------------- 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 33f3dc53e..aacd2e0cf 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -101,12 +101,10 @@ the point. The following figure represents three vectors in this manner ```{code-cell} ipython -fig, ax = plt.subplots(figsize=(10, 8)) -# Set the axes through the origin -for spine in ['left', 'bottom']: - ax.spines[spine].set_position('zero') -for spine in ['right', 'top']: - ax.spines[spine].set_color('none') +fig, ax = plt.subplots() +# Draw the axes through the origin +ax.axhline(0, color='black', lw=0.8) +ax.axvline(0, color='black', lw=0.8) ax.set(xlim=(-5, 5), ylim=(-5, 5)) ax.grid() @@ -171,12 +169,10 @@ Scalar multiplication is illustrated in the next figure --- tags: [output_scroll] --- -fig, ax = plt.subplots(figsize=(10, 8)) -# Set the axes through the origin -for spine in ['left', 'bottom']: - ax.spines[spine].set_position('zero') -for spine in ['right', 'top']: - ax.spines[spine].set_color('none') +fig, ax = plt.subplots() +# Draw the axes through the origin +ax.axhline(0, color='black', lw=0.8) +ax.axvline(0, color='black', lw=0.8) ax.set(xlim=(-5, 5), ylim=(-5, 5)) x = (2, 2) @@ -295,7 +291,7 @@ The span is a two-dimensional plane passing through these two points and the ori --- tags: [output_scroll] --- -ax = plt.figure(figsize=(10, 8)).add_subplot(projection='3d') +ax = plt.figure().add_subplot(projection='3d') x_min, x_max = -5, 5 y_min, y_max = -5, 5 @@ -692,14 +688,13 @@ x = np.linspace(xmin, xmax, 160) y = f(x) ya, yb = np.min(y), np.max(y) +# Extra height keeps the labels clear of the axes in both panels. fig, axes = plt.subplots(2, 1, figsize=(10, 10)) for ax in axes: - # Set the axes through the origin - for spine in ['left', 'bottom']: - ax.spines[spine].set_position('zero') - for spine in ['right', 'top']: - ax.spines[spine].set_color('none') + # Draw the axes through the origin + ax.axhline(0, color='black', lw=0.8) + ax.axvline(0, color='black', lw=0.8) ax.set(ylim=(-0.6, 3.2), xlim=(xmin, xmax), yticks=(), xticks=()) @@ -948,12 +943,10 @@ A = np.array(A) evals, evecs = eig(A) evecs = evecs[:, 0], evecs[:, 1] -fig, ax = plt.subplots(figsize=(10, 8)) -# Set the axes through the origin -for spine in ['left', 'bottom']: - ax.spines[spine].set_position('zero') -for spine in ['right', 'top']: - ax.spines[spine].set_color('none') +fig, ax = plt.subplots() +# Draw the axes through the origin +ax.axhline(0, color='black', lw=0.8) +ax.axvline(0, color='black', lw=0.8) ax.grid(alpha=0.4) xmin, xmax = -3, 3 From 5f9e3a873e83f68dca7750af630ae02e6b918338 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 11:59:28 +1000 Subject: [PATCH 4/9] update writing and math according to style guide --- lectures/linear_algebra.md | 154 ++++++++++++++++++++----------------- 1 file changed, 82 insertions(+), 72 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index aacd2e0cf..f1520f207 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -119,7 +119,7 @@ for v in vecs: plt.show() ``` -### Vector Operations +### Vector operations ```{index} single: Vectors; Operations ``` @@ -213,7 +213,7 @@ x + y 4 * x ``` -### Inner Product and Norm +### Inner product and norm ```{index} single: Vectors; Inner Product ``` @@ -224,7 +224,7 @@ x + y The **inner product** of vectors $x,y \in \mathbb R ^n$ is defined as $$ -x' y := \sum_{i=1}^n x_i y_i +x^\top y := \sum_{i=1}^n x_i y_i $$ Two vectors are called **orthogonal** if their inner product is zero. @@ -232,7 +232,7 @@ Two vectors are called **orthogonal** if their inner product is zero. The **norm** of a vector $x$ represents its "length" (i.e., its distance from the zero vector) and is defined as $$ -\| x \| := \sqrt{x' x} := \left( \sum_{i=1}^n x_i^2 \right)^{1/2} +\| x \| := \sqrt{x^\top x} := \left( \sum_{i=1}^n x_i^2 \right)^{1/2} $$ The expression $\| x - y\|$ is thought of as the distance between $x$ and $y$. @@ -345,7 +345,7 @@ plt.show() If $A$ contains only one vector $a_1 \in \mathbb R ^2$, then its span is just the scalar multiples of $a_1$, which is the unique line passing through both $a_1$ and the origin. -If $A = \{e_1, e_2, e_3\}$ consists of the **canonical basis vectors** of $\mathbb R ^3$, that is +If $A = \{e_1, e_2, e_3\}$ consists of the **canonical basis vectors** of $\mathbb R ^3$, that is $$ e_1 := @@ -384,7 +384,7 @@ If $y = (y_1, y_2, y_3)$ is any linear combination of these vectors, then $y_3 = Hence $A_0$ fails to span all of $\mathbb R ^3$. (la_li)= -### Linear Independence +### Linear independence ```{index} single: Vectors; Linear Independence ``` @@ -419,7 +419,7 @@ The following statements are equivalent to linear independence of $A := \{a_1, \ (The zero in the first expression is the origin of $\mathbb R ^n$) (la_unique_reps)= -### Unique Representations +### Unique representations Another nice thing about sets of linearly independent vectors is that each element in the span has a unique representation as a linear combination of these vectors. @@ -471,7 +471,7 @@ If $n = k$, then $A$ is called **square**. The matrix formed by replacing $a_{ij}$ by $a_{ji}$ for every $i$ and $j$ is called the **transpose** of $A$ and denoted $A'$ or $A^{\top}$. -If $A = A'$, then $A$ is called **symmetric**. +If $A = A^\top$, then $A$ is called **symmetric**. For a square matrix $A$, the $i$ elements of the form $a_{ii}$ for $i=1,\ldots,n$ are called the **principal diagonal**. @@ -479,7 +479,7 @@ $A$ is called **diagonal** if the only nonzero entries are on the principal diag If, in addition to being diagonal, each element along the principal diagonal is equal to 1, then $A$ is called the **identity matrix** and denoted by $I$. -### Matrix Operations +### Matrix operations ```{index} single: Matrix; Operations ``` @@ -630,7 +630,7 @@ In particular, `A @ B` is matrix multiplication, whereas `A * B` is element-by-e See [here](https://python-programming.quantecon.org/numpy.html#matrix-multiplication) for more discussion. (la_linear_map)= -### Matrices as Maps +### Matrices as maps ```{index} single: Matrix; Maps ``` @@ -649,7 +649,7 @@ You can check that this holds for the function $f(x) = A x + b$ when $b$ is the In fact, it's [known](https://en.wikipedia.org/wiki/Linear_map#Matrices) that $f$ is linear if and *only if* there exists a matrix $A$ such that $f(x) = Ax$ for all $x$. -## Solving Systems of Equations +## Solving systems of equations ```{index} single: Matrix; Solving Systems of Equations ``` @@ -747,7 +747,7 @@ A happy fact is that linear independence of the columns of $A$ also gives us uni Indeed, it follows from our {ref}`earlier discussion ` that if $\{a_1, \ldots, a_k\}$ are linearly independent and $y = Ax = x_1 a_1 + \cdots + x_k a_k$, then no $z \not= x$ satisfies $y = Az$. -### The Square Matrix Case +### The square matrix case Let's discuss some more details, starting with the case where $A$ is $n \times n$. @@ -770,7 +770,7 @@ In particular, the following are equivalent The property of having linearly independent columns is sometimes expressed as having **full column rank**. -#### Inverse Matrices +#### Inverse matrices ```{index} single: Matrix; Inverse ``` @@ -806,7 +806,7 @@ Perhaps the most important fact about determinants is that $A$ is nonsingular if This gives us a useful one-number summary of whether or not a square matrix can be inverted. -### More Rows than Columns +### More rows than columns This is the $n \times k$ case with $n > k$. @@ -839,9 +839,9 @@ $x$ that makes the distance $\| y - Ax\|$ as small as possible. To solve this problem, one can use either calculus or the theory of orthogonal projections. -The solution is known to be $\hat x = (A'A)^{-1}A'y$ --- see for example chapter 3 of [these notes](https://python.quantecon.org/_static/lecture_specific/linear_algebra/course_notes.pdf). +The solution is known to be $\hat x = (A^\top A)^{-1}A^\top y$ --- see for example chapter 3 of [these notes](https://python.quantecon.org/_static/lecture_specific/linear_algebra/course_notes.pdf). -### More Columns than Rows +### More columns than rows This is the $n \times k$ case with $n < k$, so there are fewer equations than unknowns. @@ -871,7 +871,7 @@ $$ In other words, uniqueness fails. -### Linear Equations with SciPy +### Linear equations with SciPy ```{index} single: Linear Algebra; SciPy ``` @@ -905,10 +905,10 @@ Observe how we can solve for $x = A^{-1} y$ by either via `inv(A) @ y`, or using The latter method uses a different algorithm (LU decomposition) that is numerically more stable, and hence should almost always be preferred. -To obtain the least-squares solution $\hat x = (A'A)^{-1}A'y$, use `scipy.linalg.lstsq(A, y)`. +To obtain the least-squares solution $\hat x = (A^\top A)^{-1}A^\top y$, use `scipy.linalg.lstsq(A, y)`. (la_eigen)= -## {index}`Eigenvalues ` and {index}`Eigenvectors ` +## {index}`Eigenvalues ` and {index}`eigenvectors ` ```{index} single: Linear Algebra; Eigenvalues ``` @@ -996,7 +996,7 @@ plane, although some might be repeated. Some nice facts about the eigenvalues of a square matrix $A$ are as follows -1. The determinant of $A$ equals the product of the eigenvalues. +1. The determinant of $A$ equals the product of the eigenvalues. 1. The trace of $A$ (the sum of the elements on the principal diagonal) equals the sum of the eigenvalues. 1. If $A$ is symmetric, then all of its eigenvalues are real. 1. If $A$ is invertible and $\lambda_1, \ldots, \lambda_n$ are its eigenvalues, then the eigenvalues of $A^{-1}$ are $1/\lambda_1, \ldots, 1/\lambda_n$. @@ -1025,7 +1025,7 @@ Since any scalar multiple of an eigenvector is an eigenvector with the same eigenvalue (check it), the eig routine normalizes the length of each eigenvector to one. -### Generalized Eigenvalues +### Generalized eigenvalues It is sometimes useful to consider the **generalized eigenvalue problem**, which, for given matrices $A$ and $B$, seeks generalized eigenvalues @@ -1041,12 +1041,12 @@ Of course, if $B$ is square and invertible, then we can treat the generalized eigenvalue problem as an ordinary eigenvalue problem $B^{-1} A v = \lambda v$, but this is not always the case. -## Further Topics +## Further topics We round out our discussion by briefly mentioning several other important topics. -### Series Expansions +### Series expansions ```{index} single: Linear Algebra; Series Expansions ``` @@ -1057,7 +1057,7 @@ that if $|a| < 1$, then $\sum_{k=0}^{\infty} a^k = (1 - a)^{-1}$. A generalization of this idea exists in the matrix setting. (la_mn)= -#### Matrix Norms +#### Matrix norms ```{index} single: Linear Algebra; Matrix Norms ``` @@ -1075,7 +1075,7 @@ the left-hand side is a **matrix norm** --- in this case, the so-called For example, for a square matrix $S$, the condition $\| S \| < 1$ means that $S$ is **contractive**, in the sense that it pulls all vectors towards the origin [^cfn]. (la_neumann)= -#### {index}`Neumann's Theorem ` +#### {index}`Neumann's theorem ` ```{index} single: Linear Algebra; Neumann's Theorem ``` @@ -1094,7 +1094,7 @@ $k \in \mathbb{N}$, then $I - A$ is invertible, and ``` (la_neumann_remarks)= -#### {index}`Spectral Radius ` +#### {index}`Spectral radius ` ```{index} single: Linear Algebra; Spectral Radius ``` @@ -1112,7 +1112,7 @@ there exists a $k$ with $\| A^k \| < 1$. In which case {eq}`la_neumann` is valid. -### {index}`Positive Definite Matrices ` +### {index}`Positive definite matrices ` ```{index} single: Linear Algebra; Positive Definite Matrices ``` @@ -1121,8 +1121,8 @@ Let $A$ be a symmetric $n \times n$ matrix. We say that $A$ is -1. **positive definite** if $x' A x > 0$ for every $x \in \mathbb R ^n \setminus \{0\}$ -1. **positive semi-definite** or **nonnegative definite** if $x' A x \geq 0$ for every $x \in \mathbb R ^n$ +1. **positive definite** if $x^\top A x > 0$ for every $x \in \mathbb R ^n \setminus \{0\}$ +1. **positive semi-definite** or **nonnegative definite** if $x^\top A x \geq 0$ for every $x \in \mathbb R ^n$ Analogous definitions exist for negative definite and negative semi-definite matrices. @@ -1131,12 +1131,14 @@ are strictly positive, and hence $A$ is invertible (with positive definite inverse). (la_mcalc)= -### Differentiating Linear and Quadratic Forms +### Differentiating linear and quadratic forms ```{index} single: Linear Algebra; Differentiating Linear and Quadratic Forms ``` -The following formulas are useful in many economic contexts. Let +The following formulas are useful in many economic contexts. + +Let * $z, x$ and $a$ all be $n \times 1$ vectors * $A$ be an $n \times n$ matrix @@ -1144,15 +1146,15 @@ The following formulas are useful in many economic contexts. Let Then -1. $\frac{\partial a' x}{\partial x} = a$ -1. $\frac{\partial A x}{\partial x} = A'$ -1. $\frac{\partial x'A x}{\partial x} = (A + A') x$ -1. $\frac{\partial y'B z}{\partial y} = B z$ -1. $\frac{\partial y'B z}{\partial B} = y z'$ +1. $\frac{\partial a^\top x}{\partial x} = a$ +1. $\frac{\partial A x}{\partial x} = A^\top$ +1. $\frac{\partial x^\top A x}{\partial x} = (A + A^\top) x$ +1. $\frac{\partial y^\top B z}{\partial y} = B z$ +1. $\frac{\partial y^\top B z}{\partial B} = y z^\top$ {ref}`la_ex1` below asks you to apply these formulas. -### Further Reading +### Further reading The documentation of the `scipy.linalg` submodule can be found [here](https://docs.scipy.org/doc/scipy/reference/linalg.html). @@ -1171,7 +1173,7 @@ is {cite}`Janich1994`. Let $x$ be a given $n \times 1$ vector and consider the problem $$ -v(x) = \max_{y,u} \left\{ - y'P y - u' Q u \right\} +v(x) = \max_{y,u} \left\{ - y^\top P y - u^\top Q u \right\} $$ subject to the linear constraint @@ -1191,7 +1193,7 @@ Here One way to solve the problem is to form the Lagrangian $$ -\mathcal L = - y' P y - u' Q u + \lambda' \left[A x + B u - y\right] +\mathcal L = - y^\top P y - u^\top Q u + \lambda^\top \left[A x + B u - y\right] $$ where $\lambda$ is an $n \times 1$ vector of Lagrange multipliers. @@ -1201,13 +1203,15 @@ Try applying the formulas given above for differentiating quadratic and linear f Show that these conditions imply that 1. $\lambda = - 2 P y$. -1. The optimizing choice of $u$ satisfies $u = - (Q + B' P B)^{-1} B' P A x$. -1. The function $v$ satisfies $v(x) = - x' \tilde P x$ where $\tilde P = A' P A - A'P B (Q + B'P B)^{-1} B' P A$. +1. The optimizing choice of $u$ satisfies $u = - (Q + B^\top P B)^{-1} B^\top P A x$. +1. The function $v$ satisfies $v(x) = - x^\top \tilde P x$ where $\tilde P = A^\top P A - A^\top P B (Q + B^\top P B)^{-1} B^\top P A$. As we will see, in economic contexts Lagrange multipliers often are shadow prices. ```{note} -If we don't care about the Lagrange multipliers, we can substitute the constraint into the objective function, and then just maximize $-(Ax + Bu)'P (Ax + Bu) - u' Q u$ with respect to $u$. You can verify that this leads to the same maximizer. +If we don't care about the Lagrange multipliers, we can substitute the constraint into the objective function, and then just maximize $-(Ax + Bu)^\top P (Ax + Bu) - u^\top Q u$ with respect to $u$. + +You can verify that this leads to the same maximizer. ``` ```{exercise-end} @@ -1220,7 +1224,7 @@ If we don't care about the Lagrange multipliers, we can substitute the constrain We have an optimization problem: $$ -v(x) = \max_{y,u} \{ -y'Py - u'Qu \} +v(x) = \max_{y,u} \{ -y^\top Py - u^\top Qu \} $$ s.t. @@ -1241,7 +1245,7 @@ with primitives The associated Lagrangian is: $$ -L = -y'Py - u'Qu + \lambda' \lbrack Ax + Bu - y \rbrack +L = -y^\top Py - u^\top Qu + \lambda^\top \lbrack Ax + Bu - y \rbrack $$ **Step 1.** @@ -1250,7 +1254,7 @@ Differentiating Lagrangian equation w.r.t y and setting its derivative equal to zero yields $$ -\frac{ \partial L}{\partial y} = - (P + P') y - \lambda = - 2 P y - \lambda = 0 \:, +\frac{ \partial L}{\partial y} = - (P + P^\top) y - \lambda = - 2 P y - \lambda = 0 \:, $$ since P is symmetric. @@ -1267,24 +1271,24 @@ Differentiating Lagrangian equation w.r.t. u and setting its derivative equal to zero yields $$ -\frac{ \partial L}{\partial u} = - (Q + Q') u - B'\lambda = - 2Qu + B'\lambda = 0 \: +\frac{ \partial L}{\partial u} = - (Q + Q^\top) u - B^\top\lambda = - 2Qu + B^\top\lambda = 0 \: $$ Substituting $\lambda = -2 P y$ gives $$ -Qu + B'Py = 0 \: +Qu + B^\top Py = 0 \: $$ Substituting the linear constraint $y = Ax + Bu$ into above equation gives $$ -Qu + B'P(Ax + Bu) = 0 +Qu + B^\top P(Ax + Bu) = 0 $$ $$ -(Q + B'PB)u + B'PAx = 0 +(Q + B^\top PB)u + B^\top PAx = 0 $$ which is the first-order condition for maximizing $L$ w.r.t. $u$. @@ -1292,7 +1296,7 @@ which is the first-order condition for maximizing $L$ w.r.t. $u$. Thus, the optimal choice of u must satisfy $$ -u = -(Q + B'PB)^{-1}B'PAx \:, +u = -(Q + B^\top PB)^{-1}B^\top PAx \:, $$ which follows from the definition of the first-order conditions for @@ -1304,66 +1308,66 @@ Rewriting our problem by substituting the constraint into the objective function, we get $$ -v(x) = \max_{u} \{ -(Ax+ Bu)'P(Ax+Bu) - u'Qu \} \: +v(x) = \max_{u} \{ -(Ax+ Bu)^\top P(Ax+Bu) - u^\top Qu \} \: $$ Since we know the optimal choice of u satisfies $u = -(Q + -B'PB)^{-1}B'PAx$, then +B^\top PB)^{-1}B^\top PAx$, then $$ -v(x) = -(Ax+ B u)'P(Ax+B u) - u'Q u \,\,\,\, with \,\,\,\, u = -(Q + B'PB)^{-1}B'PAx +v(x) = -(Ax+ B u)^\top P(Ax+B u) - u^\top Q u \,\,\,\, with \,\,\,\, u = -(Q + B^\top PB)^{-1}B^\top PAx $$ To evaluate the function $$ \begin{aligned} -v(x) &= -(Ax+ B u)'P(Ax+Bu) - u'Q u \\ -&= -(x'A' + u'B')P(Ax+Bu) - u'Q u \\ -&= - x'A'PAx - u'B'PAx - x'A'PBu - u'B'PBu - u'Qu \\ -&= - x'A'PAx - 2u'B'PAx - u'(Q + B'PB) u +v(x) &= -(Ax+ B u)^\top P(Ax+Bu) - u^\top Q u \\ +&= -(x^\top A^\top + u^\top B^\top)P(Ax+Bu) - u^\top Q u \\ +&= - x^\top A^\top PAx - u^\top B^\top PAx - x^\top A^\top PBu - u^\top B^\top PBu - u^\top Qu \\ +&= - x^\top A^\top PAx - 2u^\top B^\top PAx - u^\top(Q + B^\top PB) u \end{aligned} $$ -For simplicity, denote by $S := (Q + B'PB)^{-1} B'PA$, then $u = -Sx$. +For simplicity, denote by $S := (Q + B^\top PB)^{-1} B^\top PA$, then $u = -Sx$. -Regarding the second term $- 2u'B'PAx$, +Regarding the second term $- 2u^\top B^\top PAx$, $$ \begin{aligned} --2u'B'PAx &= -2 x'S'B'PAx \\ -& = 2 x'A'PB( Q + B'PB)^{-1} B'PAx +-2u^\top B^\top PAx &= -2 x^\top S^\top B^\top PAx \\ +& = 2 x^\top A^\top PB( Q + B^\top PB)^{-1} B^\top PAx \end{aligned} $$ -Notice that the term $(Q + B'PB)^{-1}$ is symmetric as both P and Q +Notice that the term $(Q + B^\top PB)^{-1}$ is symmetric as both P and Q are symmetric. -Regarding the third term $- u'(Q + B'PB) u$, +Regarding the third term $- u^\top(Q + B^\top PB) u$, $$ \begin{aligned} --u'(Q + B'PB) u &= - x'S' (Q + B'PB)Sx \\ -&= -x'A'PB(Q + B'PB)^{-1}B'PAx +-u^\top(Q + B^\top PB) u &= - x^\top S^\top (Q + B^\top PB)Sx \\ +&= -x^\top A^\top PB(Q + B^\top PB)^{-1}B^\top PAx \end{aligned} $$ Hence, the summation of second and third terms is -$x'A'PB(Q + B'PB)^{-1}B'PAx$. +$x^\top A^\top PB(Q + B^\top PB)^{-1}B^\top PAx$. This implies that $$ \begin{aligned} - v(x) &= - x'A'PAx - 2u'B'PAx - u'(Q + B'PB) u\\ - &= - x'A'PAx + x'A'PB(Q + B'PB)^{-1}B'PAx \\ - &= -x'[A'PA - A'PB(Q + B'PB)^{-1}B'PA] x + v(x) &= - x^\top A^\top PAx - 2u^\top B^\top PAx - u^\top(Q + B^\top PB) u\\ + &= - x^\top A^\top PAx + x^\top A^\top PB(Q + B^\top PB)^{-1}B^\top PAx \\ + &= -x^\top[A^\top PA - A^\top PB(Q + B^\top PB)^{-1}B^\top PA] x \end{aligned} $$ Therefore, the solution to the optimization problem -$v(x) = -x' \tilde{P}x$ follows the above result by denoting -$\tilde{P} := A'PA - A'PB(Q + B'PB)^{-1}B'PA$ +$v(x) = -x^\top \tilde{P}x$ follows the above result by denoting +$\tilde{P} := A^\top PA - A^\top PB(Q + B^\top PB)^{-1}B^\top PA$ ```{solution-end} ``` @@ -1371,4 +1375,10 @@ $\tilde{P} := A'PA - A'PB(Q + B'PB)^{-1}B'PA$ [^fn_mdt]: Although there is a specialized matrix data type defined in NumPy, it's more standard to work with ordinary NumPy arrays. See [this discussion](https://python-programming.quantecon.org/numpy.html#matrix-multiplication). -[^cfn]: Suppose that $\|S \| < 1$. Take any nonzero vector $x$, and let $r := \|x\|$. We have $\| Sx \| = r \| S (x/r) \| \leq r \| S \| < r = \| x\|$. Hence every point is pulled towards the origin. +[^cfn]: Suppose that $\|S \| < 1$. + + Take any nonzero vector $x$, and let $r := \|x\|$. + + We have $\| Sx \| = r \| S (x/r) \| \leq r \| S \| < r = \| x\|$. + + Hence every point is pulled towards the origin. From 3172e5a311f306bdbf830fe544dd7bffaa66a079 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 13:20:39 +1000 Subject: [PATCH 5/9] update the vector and matrix section --- lectures/linear_algebra.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index f1520f207..b7469c9fa 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -342,7 +342,7 @@ plt.show() #### Examples -If $A$ contains only one vector $a_1 \in \mathbb R ^2$, then its +If $A$ contains only one nonzero vector $a_1 \in \mathbb R ^2$, then its span is just the scalar multiples of $a_1$, which is the unique line passing through both $a_1$ and the origin. If $A = \{e_1, e_2, e_3\}$ consists of the **canonical basis vectors** of $\mathbb R ^3$, that is @@ -414,7 +414,7 @@ $m > n$ vectors in $\mathbb R ^n$ must be linearly dependent. The following statements are equivalent to linear independence of $A := \{a_1, \ldots, a_k\} \subset \mathbb R ^n$ 1. No vector in $A$ can be formed as a linear combination of the other elements. -1. If $\beta_1 a_1 + \cdots \beta_k a_k = 0$ for scalars $\beta_1, \ldots, \beta_k$, then $\beta_1 = \cdots = \beta_k = 0$. +1. If $\beta_1 a_1 + \cdots + \beta_k a_k = 0$ for scalars $\beta_1, \ldots, \beta_k$, then $\beta_1 = \cdots = \beta_k = 0$. (The zero in the first expression is the origin of $\mathbb R ^n$) @@ -427,13 +427,13 @@ In other words, if $A := \{a_1, \ldots, a_k\} \subset \mathbb R ^n$ is linearly independent and $$ -y = \beta_1 a_1 + \cdots \beta_k a_k +y = \beta_1 a_1 + \cdots + \beta_k a_k $$ then no other coefficient sequence $\gamma_1, \ldots, \gamma_k$ will produce the same vector $y$. -Indeed, if we also have $y = \gamma_1 a_1 + \cdots \gamma_k a_k$, +Indeed, if we also have $y = \gamma_1 a_1 + \cdots + \gamma_k a_k$, then $$ @@ -473,7 +473,7 @@ The matrix formed by replacing $a_{ij}$ by $a_{ji}$ for every $i$ and $j$ is cal If $A = A^\top$, then $A$ is called **symmetric**. -For a square matrix $A$, the $i$ elements of the form $a_{ii}$ for $i=1,\ldots,n$ are called the **principal diagonal**. +For an $n \times n$ matrix $A$, the $n$ entries $a_{11}, \ldots, a_{nn}$ form the **principal diagonal**. $A$ is called **diagonal** if the only nonzero entries are on the principal diagonal. @@ -578,7 +578,7 @@ If $I$ is the $n \times n$ identity matrix, then $IA = A$. ### Matrices in NumPy -```{index} single: Matrix; Numpy +```{index} single: Matrix; NumPy ``` NumPy arrays are also used as matrices, and have fast, efficient functions and methods for all the standard matrix operations [^fn_mdt]. From f27a03d56df3dbf1af3588fad3a600646ad14053 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Wed, 9 Sep 2026 13:32:24 +1000 Subject: [PATCH 6/9] update the eigenvector and further reading section --- lectures/linear_algebra.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index b7469c9fa..934ee83a3 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -850,7 +850,7 @@ In this case there are either no solutions or infinitely many --- in other words For example, consider the case where $k=3$ and $n=2$. -Thus, the columns of $A$ consists of 3 vectors in $\mathbb R ^2$. +Thus, the columns of $A$ consist of 3 vectors in $\mathbb R ^2$. This set can never be linearly independent, since it is possible to find two vectors that span $\mathbb R ^2$. @@ -901,9 +901,9 @@ A @ x # Should equal y solve(A, y) # Produces the same solution ``` -Observe how we can solve for $x = A^{-1} y$ by either via `inv(A) @ y`, or using `solve(A, y)`. +Observe how we can solve for $x = A^{-1} y$ either via `inv(A) @ y`, or using `solve(A, y)`. -The latter method uses a different algorithm (LU decomposition) that is numerically more stable, and hence should almost always be preferred. +To solve $Ax = y$, we generally prefer `solve(A, y)` because it avoids the unnecessary step of computing the inverse of $A$. To obtain the least-squares solution $\hat x = (A^\top A)^{-1}A^\top y$, use `scipy.linalg.lstsq(A, y)`. @@ -1021,7 +1021,7 @@ evecs Note that the *columns* of `evecs` are the eigenvectors. -Since any scalar multiple of an eigenvector is an eigenvector with the same +Since any nonzero scalar multiple of an eigenvector is an eigenvector with the same eigenvalue (check it), the eig routine normalizes the length of each eigenvector to one. @@ -1158,7 +1158,7 @@ Then The documentation of the `scipy.linalg` submodule can be found [here](https://docs.scipy.org/doc/scipy/reference/linalg.html). -Chapters 2 and 3 of the [Econometric Theory](https://johnstachurski.net/emet.html) contains +Chapters 2 and 3 of [A Primer in Econometric Theory](https://johnstachurski.net/emet) contain a discussion of linear algebra along the same lines as above, with solved exercises. If you don't mind a slightly abstract approach, a nice intermediate-level text on linear algebra From 36cf7e454523d40a49533d41027f372c5c46579a Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Thu, 10 Sep 2026 14:13:38 +1000 Subject: [PATCH 7/9] update exercise assumption and solution derivation --- lectures/linear_algebra.md | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 934ee83a3..5b26e81b8 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -1187,6 +1187,7 @@ Here * $P$ is an $n \times n$ matrix and $Q$ is an $m \times m$ matrix * $A$ is an $n \times n$ matrix and $B$ is an $n \times m$ matrix * both $P$ and $Q$ are symmetric and positive semidefinite +* $H := Q + B^\top P B$ is positive definite (What must the dimensions of $y$ and $u$ be to make this a well-posed problem?) @@ -1235,12 +1236,14 @@ $$ with primitives -- $P$ be a symmetric and positive semidefinite $n \times n$ +- $P$ is a symmetric and positive semidefinite $n \times n$ matrix -- $Q$ be a symmetric and positive semidefinite $m \times m$ +- $Q$ is a symmetric and positive semidefinite $m \times m$ matrix -- $A$ an $n \times n$ matrix -- $B$ an $n \times m$ matrix +- $A$ is an $n \times n$ matrix +- $B$ is an $n \times m$ matrix + +Assume that $H := Q + B^\top P B$ is positive definite. The associated Lagrangian is: @@ -1271,7 +1274,7 @@ Differentiating Lagrangian equation w.r.t. u and setting its derivative equal to zero yields $$ -\frac{ \partial L}{\partial u} = - (Q + Q^\top) u - B^\top\lambda = - 2Qu + B^\top\lambda = 0 \: +\frac{ \partial L}{\partial u} = - (Q + Q^\top) u + B^\top\lambda = - 2Qu + B^\top\lambda = 0 \: $$ Substituting $\lambda = -2 P y$ gives @@ -1280,7 +1283,7 @@ $$ Qu + B^\top Py = 0 \: $$ -Substituting the linear constraint $y = Ax + Bu$ into above +Substituting the linear constraint $y = Ax + Bu$ into the above equation gives $$ @@ -1315,7 +1318,7 @@ Since we know the optimal choice of u satisfies $u = -(Q + B^\top PB)^{-1}B^\top PAx$, then $$ -v(x) = -(Ax+ B u)^\top P(Ax+B u) - u^\top Q u \,\,\,\, with \,\,\,\, u = -(Q + B^\top PB)^{-1}B^\top PAx +v(x) = -(Ax+ B u)^\top P(Ax+B u) - u^\top Q u \quad \text{with} \quad u = -(Q + B^\top PB)^{-1}B^\top PAx $$ To evaluate the function @@ -1329,13 +1332,13 @@ v(x) &= -(Ax+ B u)^\top P(Ax+Bu) - u^\top Q u \\ \end{aligned} $$ -For simplicity, denote by $S := (Q + B^\top PB)^{-1} B^\top PA$, then $u = -Sx$. +For simplicity, set $S := (Q + B^\top PB)^{-1} B^\top PA$, so that $u = -Sx$. Regarding the second term $- 2u^\top B^\top PAx$, $$ \begin{aligned} --2u^\top B^\top PAx &= -2 x^\top S^\top B^\top PAx \\ +-2u^\top B^\top PAx &= +2 x^\top S^\top B^\top PAx \\ & = 2 x^\top A^\top PB( Q + B^\top PB)^{-1} B^\top PAx \end{aligned} $$ From b675a3c31d07ec4597f68b3bee9d21d1d665c769 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Fri, 11 Sep 2026 11:19:43 +1000 Subject: [PATCH 8/9] add doc link for programming series --- lectures/_config.yml | 3 +++ lectures/linear_algebra.md | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/lectures/_config.yml b/lectures/_config.yml index 92c4f0129..bd13466bf 100644 --- a/lectures/_config.yml +++ b/lectures/_config.yml @@ -104,6 +104,9 @@ sphinx: jax: - "https://jax.quantecon.org/" - null + programming: + - "https://python-programming.quantecon.org/" + - null mathjax3_config: tex: macros: diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 5b26e81b8..20373b9d0 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -65,7 +65,7 @@ These are the kinds of topics addressed by linear algebra. In this lecture we will cover the basics of linear and matrix algebra, treating both theory and computation. -We admit some overlap with [this lecture](https://python-programming.quantecon.org/numpy.html), where operations on NumPy arrays were first explained. +We admit some overlap with {doc}`programming:numpy`, where operations on NumPy arrays were first explained. Note that this lecture is more theoretical than most, and contains background material that will be used in applications as we go along. From dedc4514e976683fcfedccbefda7f38d2b0f0013 Mon Sep 17 00:00:00 2001 From: Longye Tian <133612246+longye-tian@users.noreply.github.com> Date: Fri, 11 Sep 2026 11:36:47 +1000 Subject: [PATCH 9/9] use ref link for cross-series sections --- lectures/linear_algebra.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lectures/linear_algebra.md b/lectures/linear_algebra.md index 20373b9d0..a1e3b132b 100644 --- a/lectures/linear_algebra.md +++ b/lectures/linear_algebra.md @@ -199,7 +199,7 @@ plt.show() ``` In Python, a vector can be represented as a list or tuple, such as `x = (2, 4, 6)`, but is more commonly -represented as a [NumPy array](https://python-programming.quantecon.org/numpy.html#numpy-arrays). +represented as a {ref}`NumPy array `. One advantage of NumPy arrays is that scalar multiplication and addition have very natural syntax @@ -603,13 +603,13 @@ A.shape ``` The `shape` attribute is a tuple giving the number of rows and columns --- -see [here](https://python-programming.quantecon.org/numpy.html#shape-and-dimension) +see {ref}`here ` for more discussion. To get the transpose of `A`, use `A.transpose()` or, more simply, `A.T`. There are many convenient functions for creating common matrices (matrices of zeros, -ones, etc.) --- see [here](https://python-programming.quantecon.org/numpy.html#creating-arrays). +ones, etc.) --- see {ref}`here `. Since operations are performed elementwise by default, scalar multiplication and addition have very natural syntax