06-Matrix-Multiplication.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "01777881",
   "metadata": {},
   "source": [
    "# 6 - Matrix Multiplication\n",
    "- 6.1 Standard multiplication\n",
    "- 6.2 Multiplication and equations\n",
    "- 6.3 Multiplication with diagonals\n",
    "- 6.4 LIVE EVIL\n",
    "- 6.5 Matrix-vector multiplication\n",
    "- 6.6 Creating symmetric matrices\n",
    "- 6.7 Multiply symmetric matrices\n",
    "- 6.8 Hadamard multiplication\n",
    "- 6.9 Frobenius dot product\n",
    "- 6.10 Matrix norms\n",
    "- 6.11 Matrix asymmetry index\n",
    "- 6.12 What about matrix division?\n",
    "- 6.13 Exercises\n",
    "- 6.14 Answers\n",
    "- 6.15 Code challenges\n",
    "- 6.16 Code solutions\n",
    "\n",
    "Notes, code snippets, and the end of chapter exercises from the book _Linear Algebra: Theory, Intuition, Code_ by Mike X Cohen. \n",
    "\n",
    "Find more information about the book on [github](https://github.com/mikexcohen/LinAlgBook) and [amazon](https://www.amazon.com/Linear-Algebra-Theory-Intuition-Code/dp/9083136604)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "971d2072",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebbc4e15",
   "metadata": {},
   "source": [
    "## 6.1 Standard Multiplication\n",
    "Matrix $A \\in \\mathbb{R}^{M \\times N}$ is compatible for multiplication with matrix $B \\in \\mathbb{R}^{N \\times P}$ only when the innermost dimension of $A$ and $B$ are equal.\n",
    "\n",
    "### 1/ Element Perspective\n",
    "Interpretation of matrix multiplication as dot products between rows of A and columns of B.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} & a_{1,2} & a_{1,3} \\\\\n",
    "a_{2,1} & a_{2,2} & a_{2,3} \\\\\n",
    "a_{3,1} & a_{3,2} & a_{3,3} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{1,1} & b_{1,2} \\\\\n",
    "b_{2,1} & b_{2,2} \\\\\n",
    "b_{3,1} & b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} b_{1,1} + a_{1,2} b_{2,1} + a_{1,3} b_{3,1} & a_{1,1} b_{1,2} + a_{1,2} b_{2,2} + a_{1,3} b_{3,2} \\\\\n",
    "a_{2,1} b_{1,1} + a_{2,2} b_{2,1} + a_{2,3} b_{3,1} & a_{2,1} b_{1,2} + a_{2,2} b_{2,2} + a_{2,3} b_{3,2} \\\\\n",
    "a_{3,1} b_{1,1} + a_{3,2} b_{2,1} + a_{3,3} b_{3,1} & a_{3,1} b_{1,2} + a_{3,2} b_{2,2} + a_{3,3} b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Lower triangle of product $AB$ contains dot products of rows of A whose index is larger than columns of B eg $i \\gt j$.\n",
    "- Upper triangle of product $AB$ contains dot products of columns of B whose index is larger than rows of A eg $j \\gt i$.\n",
    "- Lower and upper triangle interpretations are important for QR decomposition.\n",
    "\n",
    "### 2/ Layer Perspective\n",
    "Interpretation of matrix multiplication as sum of outer product of columns of A and rows of B.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} & a_{1,2} & a_{1,3} \\\\\n",
    "a_{2,1} & a_{2,2} & a_{2,3} \\\\\n",
    "a_{3,1} & a_{3,2} & a_{3,3} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{1,1} & b_{1,2} \\\\\n",
    "b_{2,1} & b_{2,2} \\\\\n",
    "b_{3,1} & b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} \\\\\n",
    "a_{2,1} \\\\\n",
    "a_{3,1} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{1,1} & b_{1,2} \\\\\n",
    "\\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix}\n",
    "a_{1,2} \\\\\n",
    "a_{2,2} \\\\\n",
    "a_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{2,1} & b_{2,2} \\\\\n",
    "\\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix}\n",
    "a_{1,3} \\\\\n",
    "a_{2,3} \\\\\n",
    "a_{3,3} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{3,1} & b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Layer perspective is closely related to spectral theory of matrices which says that any matrix can be represented as  sum of rank-1 matrices.\n",
    "- Spectral theory of matrices is important for singular value decomposition SVD.\n",
    "\n",
    "### 3/ Column Perspective\n",
    "Interpretation of matrix multiplication as linear weighted combination of the columns of A and the columns of B.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} & a_{1,2} & a_{1,3} \\\\\n",
    "a_{2,1} & a_{2,2} & a_{2,3} \\\\\n",
    "a_{3,1} & a_{3,2} & a_{3,3} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{1,1} & b_{1,2} \\\\\n",
    "b_{2,1} & b_{2,2} \\\\\n",
    "b_{3,1} & b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    " b_{1,1}\n",
    " \\begin{bmatrix}\n",
    " a_{1,1} \\\\\n",
    " a_{2,1} \\\\\n",
    " a_{3,1} \\\\\n",
    " \\end{bmatrix}\n",
    " +\n",
    " b_{2,1}\n",
    " \\begin{bmatrix}\n",
    " a_{1,2} \\\\\n",
    " a_{2,2} \\\\\n",
    " a_{3,2} \\\\\n",
    " \\end{bmatrix}\n",
    " +\n",
    " b_{3,1}\n",
    " \\begin{bmatrix}\n",
    " a_{1,3} \\\\\n",
    " a_{2,3} \\\\\n",
    " a_{3,3} \\\\\n",
    " \\end{bmatrix}\n",
    " &\n",
    " b_{1,2}\n",
    " \\begin{bmatrix}\n",
    " a_{1,1} \\\\\n",
    " a_{2,1} \\\\\n",
    " a_{3,1} \\\\\n",
    " \\end{bmatrix}\n",
    " +\n",
    " b_{2,2}\n",
    " \\begin{bmatrix}\n",
    " a_{1,2} \\\\\n",
    " a_{2,2} \\\\\n",
    " a_{3,2} \\\\\n",
    " \\end{bmatrix}\n",
    " +\n",
    " b_{3,2}\n",
    " \\begin{bmatrix}\n",
    " a_{1,3} \\\\\n",
    " a_{2,3} \\\\\n",
    " a_{3,3} \\\\\n",
    " \\end{bmatrix}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- The column perspective is useful in statistical infererence where the columns of the left matrix A contain the predictors and the columns of the right matrix B contain the model coefficients or weights. The product of AB is the prediction from the model.\n",
    "\n",
    "### 4/ Row Perspective\n",
    "Interpretation of matrix multiplication as linear weighted combination of the rows of A and the rows of B.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a_{1,1} & a_{1,2} & a_{1,3} \\\\\n",
    "a_{2,1} & a_{2,2} & a_{2,3} \\\\\n",
    "a_{3,1} & a_{3,2} & a_{3,3} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "b_{1,1} & b_{1,2} \\\\\n",
    "b_{2,1} & b_{2,2} \\\\\n",
    "b_{3,1} & b_{3,2} \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    " a_{1,1}\n",
    " \\begin{bmatrix}\n",
    " b_{1,1} & b_{1,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{1,2}\n",
    " \\begin{bmatrix}\n",
    " b_{2,1} & b_{2,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{1,3}\n",
    " \\begin{bmatrix}\n",
    " b_{3,1} & b_{3,2}\n",
    " \\end{bmatrix} \\\\\n",
    " a_{2,1}\n",
    " \\begin{bmatrix}\n",
    " b_{1,1} & b_{1,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{2,2}\n",
    " \\begin{bmatrix}\n",
    " b_{2,1} & b_{2,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{2,3}\n",
    " \\begin{bmatrix}\n",
    " b_{3,1} & b_{3,2}\n",
    " \\end{bmatrix} \\\\\n",
    " a_{3,1}\n",
    " \\begin{bmatrix}\n",
    " b_{1,1} & b_{1,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{3,2}\n",
    " \\begin{bmatrix}\n",
    " b_{2,1} & b_{2,2}\n",
    " \\end{bmatrix}\n",
    " +\n",
    " a_{3,3}\n",
    " \\begin{bmatrix}\n",
    " b_{3,1} & b_{3,2}\n",
    " \\end{bmatrix} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- The row perspective is useful in principal component analysis (PCA) where the rows of the right matrix B contain data (observations in rows and features in columns) and the rows of the left matrix A contain weights for combining the features.  The product of AB is the principal component scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "40fffdcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mmult_as_dot(A,B):\n",
    "    \"\"\"\n",
    "    mmult_as_dot returns the product AB as dot products of rows of A and columns of B\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: numpy.ndarray   Matrix product AB\n",
    "    \"\"\"\n",
    "    assert A.shape[1] == B.shape[0]\n",
    "    \n",
    "    m, n, p = A.shape[0], A.shape[1], B.shape[1]\n",
    "    AB = np.zeros((m,p))\n",
    "    for i in range(m):\n",
    "        for j in range(p):\n",
    "            AB[i,j] = np.dot(A[i,:], B[:,j])\n",
    "    return AB\n",
    "\n",
    "\n",
    "m, n, p = 4, 3, 2\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((n,p))\n",
    "AB = mmult_as_dot(A,B)\n",
    "expected = A @ B\n",
    "np.testing.assert_almost_equal(AB, expected, err_msg=\"mmult_as_dot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dccec95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mmult_as_outer(A,B):\n",
    "    \"\"\"\n",
    "    mmult_as_outer returns the product AB as sum of outer products of columns of A and rows of B\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: numpy.ndarray   Matrix product AB\n",
    "    \"\"\"\n",
    "    assert A.shape[1] == B.shape[0]\n",
    "    \n",
    "    m, n, p = A.shape[0], A.shape[1], B.shape[1]\n",
    "    AB = np.zeros((m,p))\n",
    "    for k in range(n):\n",
    "        AB += np.outer(A[:,k], B[k,:])\n",
    "    return AB\n",
    "\n",
    "\n",
    "m, n, p = 4, 3, 2\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((n,p))\n",
    "AB = mmult_as_outer(A,B)\n",
    "expected = A @ B\n",
    "np.testing.assert_almost_equal(AB, expected, err_msg=\"mmult_as_outer\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6bf5fd2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mmult_as_col(A,B):\n",
    "    \"\"\"\n",
    "    mmult_as_col returns the product AB as linear weighted combination of columns of A and B\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: numpy.ndarray   Matrix product AB\n",
    "    \"\"\"\n",
    "    assert A.shape[1] == B.shape[0]\n",
    "    \n",
    "    m, n, p = A.shape[0], A.shape[1], B.shape[1]\n",
    "    AB = np.empty((m,p))\n",
    "    for k in range(p):\n",
    "        AB[:,k] = A @ B[:,k]\n",
    "    return AB\n",
    "\n",
    "\n",
    "m, n, p = 4, 3, 2\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((n,p))\n",
    "AB = mmult_as_col(A,B)\n",
    "expected = A @ B\n",
    "np.testing.assert_almost_equal(AB, expected, err_msg=\"mmult_as_col\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6d4734e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mmult_as_row(A,B):\n",
    "    \"\"\"\n",
    "    mmult_as_row returns the product AB as linear weighted combination of rows of A and B\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: numpy.ndarray   Matrix product AB\n",
    "    \"\"\"\n",
    "    assert A.shape[1] == B.shape[0]\n",
    "    \n",
    "    m, n, p = A.shape[0], A.shape[1], B.shape[1]\n",
    "    AB = np.empty((m,p))\n",
    "    for i in range(m):\n",
    "        AB[i,:] = A[i,:] @ B\n",
    "    return AB\n",
    "\n",
    "\n",
    "m, n, p = 4, 3, 2\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((n,p))\n",
    "AB = mmult_as_row(A,B)\n",
    "expected = A @ B\n",
    "np.testing.assert_almost_equal(AB, expected, err_msg=\"mmult_as_row\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0e8a4c7",
   "metadata": {},
   "source": [
    "## 6.2 Multiplication and equations\n",
    "When performing matrix algebra to both sides of an equation, then matrix order eg pre-multiplication or post-multiplication must be preserved.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A &= A \\\\\n",
    "AB &= AB \\quad \\text{post-multiply} \\\\\n",
    "BA &= BA \\quad \\text{pre-multiply} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7303987",
   "metadata": {},
   "source": [
    "## 6.3 Multiplication with diagonals\n",
    "Pre-multiplication with a diagonal matrix scales the rows of the right matrix by diagonal elements.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "k_1 & 0 & 0 \\\\\n",
    "0 & k_2 & 0 \\\\\n",
    "0 & 0 & k_3 \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "g & h & i \\\\\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "k_1 a & k_1 b & k_1 c \\\\\n",
    "k_2 d & k_2 e & k_2 f \\\\\n",
    "k_3 g & k_3 h & k_3 i \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Post-multiplication with a diagonal matrix scales the columns of the left matrix by diagonal elements.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "g & h & i \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "k_1 & 0 & 0 \\\\\n",
    "0 & k_2 & 0 \\\\\n",
    "0 & 0 & k_3 \\\\\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "k_1 a & k_2 b & k_3 c \\\\\n",
    "k_1 d & k_2 e & k_3 f \\\\\n",
    "k_1 g & k_2 h & k_3 i \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Multiplication of two diagonal matrices is the product of the diagonals.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & 0 & 0 \\\\\n",
    "0 & b & 0 \\\\\n",
    "0 & 0 & c \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "d & 0 & 0 \\\\\n",
    "0 & e & 0 \\\\\n",
    "0 & 0 & f \\\\\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "a d & 0 & 0 \\\\\n",
    "0 & b e & 0 \\\\\n",
    "0 & 0 & c f \\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3ea554a",
   "metadata": {},
   "source": [
    "## 6.4 LIVE EVIL\n",
    "An operation applied to the product of matrices is equal to the product of the operation applied to each matrix but in reverse order.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "(A_1 A_2 \\cdots A_N)^T &= A_N^T \\cdots A_2^T A_1^T \\\\\n",
    "(A_1 A_2 \\cdots A_N)^{-1} &= A_N^{-1} \\cdots A_2^{-1} A_1^{-1} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d8c770",
   "metadata": {},
   "source": [
    "## 6.5 Matrix-vector multiplication\n",
    "Matrix-vector multiplication always produces a vector.\n",
    "- This is the connection between linear transformations that underpin computer graphics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08a9e2d7",
   "metadata": {},
   "source": [
    "## 6.6 Creating symmetric matrices\n",
    "Methods for creating symmetric matrices from a non-symmetric matrix.\n",
    "\n",
    "Notes\n",
    "- Symmetric matrix implies $C = C^T$.\n",
    "- In general, the matrix produced from Method 1 will not be equal to matrix produced from Method 2 eg not unique.\n",
    "\n",
    "### Method 1/ Additive\n",
    "$C = \\frac{1}{2} (A^T + A)$\n",
    "\n",
    "Notes\n",
    "- A must be square\n",
    "\n",
    "### Method 2/ Multiplicative\n",
    "$C = A^T A$\n",
    "\n",
    "Notes\n",
    "- C will be square, even if A is not"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa9e8809",
   "metadata": {},
   "source": [
    "## 6.8 Hadamard multiplication\n",
    "Hadamard multiplication is the elementwise multiplication of 2 matrices. Both matrics must be of same dimension.\n",
    "\n",
    "$$\n",
    "C = A \\odot B \\implies c_{i,j} = a_{i,j} b_{i,j}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Unlike matrix multiplication, Hadamard multiplication is commutative $A \\odot B = B \\odot A$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc8353e0",
   "metadata": {},
   "source": [
    "## 6.9 Frobenius dot product\n",
    "Frobenius dot product is the sum of the elementwise multiplication of 2 matrices. Both matrices must be of same dimension.\n",
    "\n",
    "$$\n",
    "\\langle A, B \\rangle_{F} = \\sum_i^M \\sum_j^N a_{i,j} b_{i,j} = \\sum A \\odot B\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Frobenius dot product is equivalent to the sum of the elements of the Hadamard matrix.\n",
    "* Frobenius dot product is used as a measure of distance or similarity between two matrices.\n",
    "\n",
    "### Alternate Definition\n",
    "$$\n",
    "\\langle A, B \\rangle_{F} = \\text{tr}(A^T B)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0980f7e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def frob(A,B):\n",
    "    \"\"\"\n",
    "    frob returns the frobenius dot product between matrices A and B\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: float           Frobenius dot product\n",
    "    \"\"\"\n",
    "    assert A.shape == B.shape\n",
    "    \n",
    "    return np.sum(A * B)\n",
    "\n",
    "\n",
    "m, n = 4, 3\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((m,n))\n",
    "ABF = frob(A,B)\n",
    "expected = np.trace(A.T @ B)\n",
    "np.testing.assert_almost_equal(ABF, expected, err_msg=\"frob\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3def2f13",
   "metadata": {},
   "source": [
    "## 6.10 Matrix norms\n",
    "There exists many different kinds of matrix norms.  Some of them are listed below.\n",
    "\n",
    "Frobenius matrix norm aka _L2 norm_.\n",
    "\n",
    "$$\n",
    "|| A ||_F = \\sqrt{ \\sum_i^M \\sum_j^N a_{i,j}^2 }\n",
    "$$\n",
    "\n",
    "Properties\n",
    "* Frobenius matrix norm can be used to measure distance between matrices $|| A - B ||_F$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7a804b28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def distL2(A,B):\n",
    "    \"\"\"\n",
    "    distL2 returns distance between two matrices A and B using the L2 norm\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :param B: numpy.ndarray  Matrix B\n",
    "    :return: float           Frobenius dot product\n",
    "    \"\"\"\n",
    "    assert A.shape == B.shape\n",
    "    \n",
    "    return math.sqrt(np.sum(np.square(A-B)))\n",
    "\n",
    "\n",
    "m, n = 4, 3\n",
    "A = np.random.random((m,n))\n",
    "B = np.random.random((m,n))\n",
    "dAB = distL2(A,B)\n",
    "expected = np.linalg.norm(A-B, 'fro')\n",
    "np.testing.assert_almost_equal(dAB, expected, err_msg=\"distL2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b56dc2a7",
   "metadata": {},
   "source": [
    "## 6.11 Matrix asymmetry index\n",
    "Every square matrix $A$ can be expressed as the sum of a symmetric matrix $A_s$ and skew-symmetric matrix $A_k$.\n",
    "\n",
    "$$\n",
    "A = A_s + A_k\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Skew-symmetric implies $C = -C^T$ (entries along diagonals must be 0).\n",
    "- Think of the skew-symmetric matrix $A_k$ as the residual. When this matrix is small, then the original matrix is relatively symmetric.\n",
    "\n",
    "### Skew-Symmetric Matrix, $A_k$\n",
    "A skew-symmetric matrix is formed by re-arranging the terms in the additive method for constructing a symmetric matrix.\n",
    "\n",
    "$A_k = (A - A^T) / 2$\n",
    "\n",
    "### Symmetric Matrix, $A_s$\n",
    "\n",
    "$A_s = A - A_k$\n",
    "\n",
    "### Matrix Asymmetry Index, $A_\\sigma$\n",
    "Matrix asymetry index describes how symmetric a matrix is.\n",
    "- Values near 0 indicate the matrix is (nearly) symmetric.\n",
    "- Values near 1 indicate the matrix is (nearly) skew-symmetric.\n",
    "\n",
    "$$\n",
    "A_\\sigma = \\frac{||A_k||_F^2}{||A||_F^2}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- Invert the value of index to obtain a symmetry index $1 - A_\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8ca0b047",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4\n",
    "A = np.random.random((m,m))\n",
    "\n",
    "Asym = 0.5*(A.T + A)\n",
    "Assym = 0.5*(A - A.T)\n",
    "Asum = Asym + Assym\n",
    "np.testing.assert_almost_equal(Asum, A, err_msg=\"matrix as sum of symmetric and skew-symmetric\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9f6f001",
   "metadata": {},
   "source": [
    "## 6.15 Code Challenges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47b6d09",
   "metadata": {},
   "source": [
    "> Implement matrix multiplication between a $2 \\times 4$ matrix and a $4 \\times 3$ matrix using the layer perspective."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9851fad2",
   "metadata": {},
   "source": [
    "See section 6.1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c70d7b65",
   "metadata": {},
   "source": [
    "> Generate a 4 x 4 diagonal matrix and a 4 x 4 dense matrix of random numbers. Compute both standard and Hadamard multiplication between them. You already know that the resulting product matrices are not the same, but what about the diagonals of those two product matrices?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "449df5b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AB\n",
      " [[0.34413881 0.05851533 0.25849147 0.76693329]\n",
      " [0.31091244 0.88084554 0.93684188 0.19550658]\n",
      " [0.23027404 0.44177244 1.30454044 2.67610962]\n",
      " [3.53445122 0.54135024 1.28286597 2.36713607]]\n",
      "ABpt\n",
      " [[0.34413881 0.         0.         0.        ]\n",
      " [0.         0.88084554 0.         0.        ]\n",
      " [0.         0.         1.30454044 0.        ]\n",
      " [0.         0.         0.         2.36713607]]\n"
     ]
    }
   ],
   "source": [
    "A = np.diag(np.arange(1, 5, dtype=float))\n",
    "B = np.random.random((4,4))\n",
    "AB = A @ B\n",
    "ABpt = A * B\n",
    "\n",
    "np.testing.assert_almost_equal(np.diag(AB), np.diag(ABpt), \n",
    "                               err_msg=\"diagonals of matrix and pointwise product are equal\")\n",
    "print(\"AB\\n\", AB)\n",
    "print(\"ABpt\\n\", ABpt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3851d424",
   "metadata": {},
   "source": [
    "> Consider $C_1 = (A^T + A) / 2$ and $C_2 = A^T A$ for some square non-symmetric matrix A. $C_1$ and $C_2$ are both symmetric and both formed from the same matrix, yet in general $C_1 \\neq C_2$. Interestingly, if $A$ is a diagonal matrix with all non-negative values, then $C_1 = C_2^{\\frac{1}{2}}$. Show this in code using random numbers, and then explain why you get the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1addc3fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sym_add(A):\n",
    "    \"\"\"\n",
    "    sym_add returns a symmetric matrix from A using the additive method\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :return: numpy.ndarray   Symmetric matrix A = A^T\n",
    "    \"\"\"\n",
    "    return (A + A.T) / 2.\n",
    "\n",
    "\n",
    "def sym_mult(A):\n",
    "    \"\"\"\n",
    "    sym_mult returns a symmetric matrix from A using the multiplicative method\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :return: numpy.ndarray   Symmetric matrix A = A^T\n",
    "    \"\"\"\n",
    "    return A @ A.T\n",
    "\n",
    "\n",
    "# Since A is a diagonal matrix, A @ A.T contains the square of the diagonals.\n",
    "A = np.diag(np.random.random(4))\n",
    "C1 = sym_add(A)\n",
    "C2 = sym_mult(A)\n",
    "np.testing.assert_almost_equal(C1, np.sqrt(C2), err_msg=\"C_1 = C_2^(1/2)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9016a7e",
   "metadata": {},
   "source": [
    "> Let's explore the Cauchy-Schwartz inequality. Generate a random matrix $A$ and a random vector $v$. Show that $||Av|| \\leq ||A||_F ||v||$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9635a9d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.random((4,4))\n",
    "v = np.random.random(4)\n",
    "\n",
    "Avnorm = np.linalg.norm(A @ v)\n",
    "Anorm = np.linalg.norm(A, 'fro')\n",
    "vnorm = np.linalg.norm(v)\n",
    "CS = bool(Avnorm <= Anorm * vnorm)\n",
    "np.testing.assert_equal(CS, True, err_msg=\"Cauchy-Schwartz inequality\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72ff95e4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}