`np.sum(a * b)` versus `a @ b`

[oyamad]

Inner dots are computed as np.sum(a * b) in many places. Consider replacing it with a @ b. (See also the 3rd item in #169 (comment).)

Below is the output of ChatGPT:

===

• Under the hood: a @ b (or np.dot(a, b)) delegates to optimized BLAS libraries that implement fused multiply-add (FMA) operations, eliminating the temporary array and reducing rounding errors compared to separate multiply and sum steps. (https://stackoverflow.com/questions/18092984/numpy-difference-between-dota-b-and-ab-sum)
• Performance overhead: np.sum(a * b) allocates an intermediate array and then sums it, incurring extra memory traffic and Python-level looping; by contrast, BLAS kernels stream data efficiently in optimized C code. (https://stackoverflow.com/questions/75556221/why-is-np-dot-so-much-faster-than-np-sum)
• Numerical accuracy: BLAS libraries use hardware FMA instructions to fuse multiply and add into a single rounding operation, improving precision over two separate operations. (https://stackoverflow.com/questions/45200278/numpy-fusing-multiply-and-add-to-avoid-wasting-memory)
• Official recommendation: NumPy’s documentation advises using np.dot(a, b) or the @ operator (via np.matmul, per PEP 465) for dot products and matrix multiplication. (https://numpy.org/doc/stable/reference/generated/numpy.dot.html, https://numpy.org/doc/stable/reference/generated/numpy.matmul.html, https://peps.python.org/pep-0465/)

===

[jstac]

Thanks @oyamad , this is very helpful -- and interesting.

@mmcky -- would you mind to supervise fixes? Perhaps @oyamad 's students or one of our RAs can take this job?

[mmcky]

These are the following instances using a regex pattern

• lectures/aiyagari.md: return np.sum(asset_probs * am.a_vals) [aiyagari] MAINT: transfer np.sum(a * b) to a @ b #467
• lectures/career.md: self.F_mean = np.sum(self.θ * self.F_probs)
• lectures/career.md: self.G_mean = np.sum(self.ϵ * self.G_probs)
• lectures/lake_model.md: β * γ * np.sum(np.maximum(U, V) * p_vec)
• lectures/lake_model.md: w = np.sum(V * p_vec * (w_vec - τ > w_bar)) / np.sum(p_vec * (w_vec -
• lectures/linear_algebra.md:np.sum(x * y) # Inner product of x and y
• lectures/mccall_model_with_separation.md: d_new = np.sum(np.maximum(v, u(c) + β * d) * q)
• lectures/mccall_model.md: reject = c + β * np.sum(v * q)
• lectures/mccall_model.md: return (1 - β) * (c + β * np.sum(v * q))
• lectures/mccall_model.md: h = np.sum(w * q) / (1 - β)
• lectures/mccall_model.md: h_next = c + β * np.sum(s * q)
• lectures/mccall_q.md: reject = c + β * np.sum(v * q)
• lectures/mle.md: return np.sum(y * np.log(μ) - μ - np.log(factorial(y)))
• lectures/mle.md: return np.sum(y * np.log(μ) + (1 - y) * np.log(1 - μ))
• lectures/mle.md: return np.sum((X.T * y * ϕ / μ - X.T * (1 - y) * ϕ / (1 - μ)),
• lectures/multivariate_normal.md:Σθ_hat_arr_C = np.array([np.sum(C[n, i+1:n+1] ** 2) for i in range(n)])
• lectures/opt_transport.md:total_cost = np.sum(X * C)
• lectures/prob_matrix.md:mean = np.sum(Benford_pmf * k)
• lectures/prob_matrix.md:var = np.sum([(k-mean)**2 * Benford_pmf])

[mmcky]

Suggested approach:

1. change the current code-cell to hide-input and add a cell to save the result for comparison with new code
2. add a new code-cell with the updated code
3. include another code-cell with an assert to test equality of the former result and the new result
4. ask co-pilot to provide a review (optional)

Once this has passed testing

4. Ask for review of the PR, and/or remove the original and comparison code

Example: #467

[oyamad]

On lectures/linear_algebra.md:np.sum(x * y) # Inner product of x and y:
This requires some discussion: Does the lecture want to teach element-wise multiplication (x * y) and power (x ** 2)?

If not, np.sum(x * y) should be replaced with x @ y, and np.sqrt(np.sum(x**2)) with np.sqrt(x @ x).

[oyamad]

Or maybe have both "take 1: np.sum(x * y)" and "take 2: x @ y" and explain that the latter is better?

[mmcky]

Or maybe have both "take 1: np.sum(x * y)" and "take 2: x @ y" and explain that the latter is better?

@oyamad thank you. I think this suggestion is great. It is a confusing aspect of numpy and demonstrating the difference between both seems like a good addition to me.

cc: @jstac

[mmcky]

@suda-yuga and @UdohNakamura -- thank you for your PR's ❤ -- I have enabled all the automated checks to be run and will organize to review them.