{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "###
San Jose State University
Department of Applied Data Science
\n", "#
DATA 220
Mathematical Methods for Data Analysis
\n", "###
Spring 2021
Instructor: Ron Mak
\n", "#
Assignment #11
Matrix Operations Problem Set
SOLUTION
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy.linalg import solve, inv, qr, norm\n", "\n", "%precision 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### PROBLEM 1. Amalgamated Widgets, Inc. manufactures several products using various amounts of materials per product:\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ProductMetalWoodPlasticClothPaper
A01.30.20.80.4
B001.50.40.3
C0.25000.20.7
D000.30.70.5
E0.1500.50.40.8
\n", "\n", "#### On a certain day, the following total amounts of materials were used:\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
MetalWoodPlasticClothPaper
51.0312.0215.4373.1356.0
\n", "\n", "#### How many of each product were made that day?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### SOLUTION. Let xA, xB, xC, xD, xE be the numbers made of each product. Then we must have\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0xA+0xB+0.25xC+0xD+0.15xE=51.0
1.3xA+0xB+0xC+0xD+0xE=312.0
0.2xA+1.5xB+0xC+0.3xD+0.5xE=215.4
0.8xA+0.4xB+0.2xC+0.7xD+0.4xE=373.1
0.4xA+0.3xB+0.7xC+0.5xD+0.8 xE=356.0
\n", "\n", "#### So we really want the transpose of the above materials table:\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", "
MaterialABCDE
Metal000.2500.15
Wood1.30000
Plastic0.21.500.30.5
Cloth0.80.40.20.70.4
Paper0.40.30.70.50.8
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = ['A', 'B', 'C', 'D', 'E']\n", "m = ['metal', 'wood', 'plastic', 'cloth', 'paper']\n", "\n", "# Matrix of materials used.\n", "M = np.array([[0, 1.3, 0.2, 0.8, 0.4],\n", " [0, 0, 1.5, 0.4, 0.3],\n", " [0.25, 0, 0, 0.2, 0.7],\n", " [0, 0, 0.3, 0.7, 0.5],\n", " [0.15, 0, 0.5, 0.4, 0.8]])\n", "\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "product_materials = []\n", "\n", "print('Product Metal Wood Plastic Cloth Paper')\n", "print('---------------------------------------------------')\n", "\n", "for i in range(len(p)):\n", " product_materials.append(M[i])\n", " print(f'{p[i]} ', end='')\n", " \n", " for j in range(len(m)):\n", " print(f'{product_materials[i][j]:9.2f}', end='')\n", "\n", " print()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Vector b of total amounts.\n", "b = np.array([51.0, 312.0, 215.4, 373.1, 356.0])\n", "b" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(' Metal Wood Plastic Cloth Paper')\n", "print('---------------------------------------------')\n", "\n", "for i in range(len(b)):\n", " print(f'{b[i]:9.1f}', end='')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Transpose the materials matrix.\n", "A = M.T\n", "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Vector n of numbers of products made.\n", "n = solve(A, b)\n", "n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('Product Number')\n", "print('------- ------')\n", "\n", "for i in range(len(p)):\n", " print(f'{p[i]} {n[i]:6.2f}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verification: Should equal b.\n", "A@n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### PROBLEM 2. Find the least-squares regression line for the points (1,1), (2,2), (3,4), (4,4), and (5,6) using the matrix form of the normal equations A.T@A@x = A.T@b.\n", "\n", "#### SOLUTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = np.array([[1, 1],\n", " [1, 2],\n", " [1, 3],\n", " [1, 4],\n", " [1, 5]])\n", "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b = np.array([1, 2, 4, 4, 6])\n", "b" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "AT = A.T\n", "AT" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATA = AT@A\n", "ATA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%precision 3\n", "\n", "ATAinv = inv(ATA)\n", "ATAinv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATb = AT@b\n", "ATb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Coefficients of the least-squares regression line:\")\n", "\n", "xhat = ATAinv@ATb\n", "xhat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Test the least-squares regression line: Should approximate b.\n", "# y = a0 + a1*x\n", "A@xhat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### PROBLEM 3. The following table shows the mean distances x in astonomical units (multiples of the distance between the sun and the earth) and the periods y in years (how many earth-years to orbit the sun):\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PlanetMercuryVenusEarthMarsJupiterSaturn
Distance x0.3870.7231.0001.5245.2039.537
Years y0.2410.6151.0001.88111.86229.457
\n", "\n", "#### It turns out that if you take the natural logarithm of the numbers in the table, the transformed points lie roughly in a straight line. Find the least-squares regression line using the matrix form of the normal equations A.T@A@x = A.T@b. \n", "\n", "#### SOLUTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Matrix of planetary data.\n", "P = np.array([[0.387, 0.723, 1.000, 1.524, 5.203, 9.537],\n", " [0.241, 0.615, 1.000, 1.881, 11.862, 29.457]])\n", "P" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Matrix of natural logs of planetary data.\n", "L = np.log(P)\n", "L" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Initialize matrix A using the first row of matrix L.\n", "A = L[0].reshape(len(L[0]), 1)\n", "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Insert the column of ones.\n", "A = np.insert(A, 0, 1, axis=1)\n", "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Vector b is the second row of matrix L.\n", "b = L[1]\n", "b" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "AT = A.T\n", "AT" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATA = AT@A\n", "ATA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATAinv = inv(ATA)\n", "ATAinv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATb = AT@b\n", "ATb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Coefficients of the least-squares regression line:\")\n", "\n", "xhat = ATAinv@ATb\n", "xhat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Therefore, ln y = 0 + (3/2)ln x, and y = x3/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### PROBLEM 4.a. Here's more data about that certain day at Amalgamated Widgets, Inc. There are three more products and the cost to make one of each product.\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ProductMetalWoodPlasticClothPaperCost
A01.30.20.80.421.60
B001.50.40.317.20
C0.25000.20.74.70
D000.30.70.59.60
E0.1500.50.40.815.70
F1.000.900.000.650.0026.80
G0.220.001.801.001.2525.50
H0.200.750.651.001.2522.40
\n", "\n", "#### The manufacturing cost of a product is related to the amounts of materials used to make the product. Find the least-squares regression line using the matrix form of the normal equations A.T@A@x = A.T@b that will compute an estimated manufacturing cost of a product.\n", "\n", "#### SOLUTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']\n", "m = ['metal', 'wood', 'plastic', 'cloth', 'paper']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Matrix of materials used.\n", "M = np.array([[0, 1.3, 0.2, 0.8, 0.4],\n", " [0, 0, 1.5, 0.4, 0.3],\n", " [0.25, 0, 0, 0.2, 0.7],\n", " [0, 0, 0.3, 0.7, 0.5],\n", " [0.15, 0, 0.5, 0.4, 0.8],\n", " [1.0, 0.9, 0, 0.65, 0],\n", " [0.22, 0, 1.8, 1.0, 1.25],\n", " [0.2, 0.75, 0.65, 1.0, 1.25]])\n", "\n", "M" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "product_materials = []\n", "\n", "print('Product Metal Wood Plastic Cloth Paper')\n", "print('---------------------------------------------------')\n", "\n", "for i in range(len(p)):\n", " product_materials.append(M[i])\n", " print(f'{p[i]} ', end='')\n", " \n", " for j in range(len(m)):\n", " print(f'{product_materials[i][j]:9.2f}', end='')\n", "\n", " print()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Cost vector b.\n", "b = np.array([21.60, 17.20, 4.70, 9.60, 15.70, 26.80, 25.50, 22.40])\n", "\n", "print('Product Cost')\n", "print('--------------')\n", "\n", "for i in range(len(p)):\n", " print(f'{p[i]} {b[i]:7.2f}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Matrix A is matrix M with an inserted column of ones.\n", "A = np.insert(M, 0, 1, axis=1)\n", "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "AT = A.T\n", "AT" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATA = AT@A\n", "ATA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATAinv = inv(ATA)\n", "ATAinv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ATb = AT@b\n", "ATb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Coefficients of the least-squares regression line:\")\n", "\n", "xhat = ATAinv@ATb\n", "xhat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Test the least-squares regression line: Should approximate b.\n", "A@xhat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "#### PROBLEM 4.b. Find the least-squares regression line for the same data using QR factorization. Verify that matrix R is upper-triangular and invertible. Verify that the columns of matrix Q are orthonormal. \n", "\n", "#### SOLUTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Q,R = qr(A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verify R is upper-triangular.\n", "R" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verify R is invertible.\n", "Rinv = inv(R)\n", "Rinv" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verify R is invertible: Should be the identify matrix I.\n", "Rinv@R" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Q" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "QT = Q.T\n", "QT" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verify Q is orthonormal: The norms of each column of Q \n", "# (here using the rows of Q.T) should be 1.\n", "for r in QT:\n", " print(norm(r))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Verify Q is orthonormal: The dot product of each column of Q \n", "# (here using the rows of Q.T) \n", "# with each other column should be 0.\n", "for i in range(len(QT) - 1):\n", " for j in range(i + 1, len(QT)):\n", " print(f'columns {i} and {j}: {QT[i]@QT[j]}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Should be the identify matrix I.\n", "QT@Q" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Coefficients of the least-squares regression line:\")\n", "\n", "xhat = Rinv@QT@b\n", "xhat" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }