- "markdown": "# Matrix Calculations {#sec-matrix-calcs}\n\n::: callout-advanced\nThis entire section is only appropriate if you've already had linear algebra.\n:::\n\n## {{< fa bullseye >}} Objectives\n\n- Understand how to do matrix algebra in relevant programming languages\n\n\n\nWhile R, SAS, and Python are all extremely powerful statistical programming languages, the core of most programming languages is the ability to do basic calculations and matrix arithmetic. \nAs almost every dataset is stored as a matrix-like structure (data sets and data frames both allow for multiple types, which isn't quite compatible with more canonical matrices), it is useful to know how to do matrix-level calculations in whatever language you are planning to use to work with data.\n\nIn this section, we will essentially be using our programming language as overgrown calculators.\n\n## Matrix Operations\n\n| Operation | R | SAS | Python |\n|:---------|:-------|:-------|:-------|\n| Addition | \\+ | \\+ | \\+ |\n| Subtraction | \\- | \\- | \\- |\n| Elementwise Multiplication | \\* | \\# | \\* |\n| Division | / | / | / |\n| Modulo (Remainder) | %% | MOD | % |\n| Integer Division | %/% | FLOOR(x\\\\y) | // |\n| Elementwise Exponentiation | \\^ | \\## | `**` |\n| Matrix/Vector Multiplication | %\\*% | \\* | `np.dot()` |\n| Matrix Exponentiation | \\^ | \\*\\* | `np.exp()` |\n| Matrix Transpose | `t(A)` | `A`\\` | `np.transpose(A)` | |\n| Matrix Determinant | `det(A)` | `det(A)` | `np.linalg.det(A)` |\n| Matrix Diagonal | `diag(A)` | `diag(A)` | `np.linalg.diag(A)` |\n| Matrix Inverse | `solve(A)` | `solve(A, diag({...}))` | `np.linalg.inv(A)` |\n\n: Table of common mathematical and matrix operations in R, SAS, and Python [@PythonProgramming]. {#tbl-math-ops}\n\n::: callout-demo\n### Basic Mathematical Operators \n\n::: panel-tabset\n\n#### R {.unnumbered}\n\n\n::: {.cell}\n\n```{.r .cell-code}\nx <- 1:10\ny <- seq(3, 30, by = 3)\n\nx + y\n## [1] 4 8 12 16 20 24 28 32 36 40\nx - y\n## [1] -2 -4 -6 -8 -10 -12 -14 -16 -18 -20\nx * y\n## [1] 3 12 27 48 75 108 147 192 243 300\nx / y\n## [1] 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333\n## [8] 0.3333333 0.3333333 0.3333333\nx^2\n## [1] 1 4 9 16 25 36 49 64 81 100\nt(x) %*% y\n## [,1]\n## [1,] 1155\n```\n:::\n\n\n#### Python {.unnumbered}\n\n\n::: {.cell}\n\n```{.python .cell-code}\nimport numpy as np\n\nx = np.array(range(1, 11))\ny = np.array(range(3, 33, 3)) # python indexes are not inclusive\n\nx + y\n## array([ 4, 8, 12, 16, 20, 24, 28, 32, 36, 40])\nx - y\n## array([ -2, -4, -6, -8, -10, -12, -14, -16, -18, -20])\nx * y\n## array([ 3, 12, 27, 48, 75, 108, 147, 192, 243, 300])\nx / y\n## array([0.33333333, 0.33333333, 0.33333333, 0.33333333, 0.33333333,\n## 0.33333333, 0.33333333, 0.33333333, 0.33333333, 0.33333333])\nx ** 2\n## array([ 1, 4, 9, 16, 25, 36, 49, 64, 81, 100])\nnp.dot(x.T, y)\n## 1155\n```\n:::\n\n\n#### SAS {.unnumbered}\n\nBy default, SAS creates row vectors with `do(a, b, by = c)` syntax. The transpose operator (a single backtick) can be used to transform `A` into `A`\\`.\n\n proc iml; \n x = do(1, 10, 1);\n y = do(3, 30, 3);\n\n z = x + y;\n z2 = x - y;\n z3 = x # y;\n z4 = x/y;\n z5 = x##2;\n z6 = x` * y;\n print z, z2, z3, z4, z5, z6;\n quit;\n:::\n:::\n\n::: callout-demo\n### Matrix Operations\n\nOther matrix operations, such as determinants and extraction of the matrix diagonal, are similarly easy:\n\n::: panel-tabset\n\n#### R {.unnumbered}\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmat <- matrix(c(1, 2, 3, 6, 4, 5, 7, 8, 9), nrow = 3, byrow = T)\nmat\n## [,1] [,2] [,3]\n## [1,] 1 2 3\n## [2,] 6 4 5\n## [3,] 7 8 9\nt(mat) # transpose\n## [,1] [,2] [,3]\n## [1,] 1 6 7\n## [2,] 2 4 8\n## [3,] 3 5 9\ndet(mat) # get the determinant\n## [1] 18\ndiag(mat) # get the diagonal\n## [1] 1 4 9\ndiag(diag(mat)) # get a square matrix with off-diag 0s\n## [,1] [,2] [,3]\n## [1,] 1 0 0\n## [2,] 0 4 0\n## [3,] 0 0 9\ndiag(1:3) # diag() also will create a diagonal matrix if given a vector\n## [,1] [,2] [,3]\n## [1,] 1 0 0\n## [2,] 0 2 0\n## [3,] 0 0 3\n```\n:::\n\n\n#### Python {.unnumbered}\n\n\n::: {.cell}\n\n```{.python .cell-code}\nimport numpy as np\nmat = np.array([[1, 2, 3],[6, 4, 5],[7, 8, 9]], dtype = int, order ='C')\n\nmat\n## array([[1, 2, 3],\n## [6, 4, 5],\n## [7, 8, 9]])\nmat.T\n## array([[1, 6, 7],\n## [2, 4, 8],\n## [3, 5, 9]])\nnp.linalg.det(mat) # numerical precision...\n## 18.000000000000004\nnp.diag(mat)\n## array([1, 4, 9])\nnp.diag(np.diag(mat))\n## array([[1, 0, 0],\n## [0, 4, 0],\n## [0, 0, 9]])\nnp.diag(range(1, 4))\n## array([[1, 0, 0],\n## [0, 2, 0],\n## [0, 0, 3]])\n```\n:::\n\n\n#### SAS {.unnumbered}\n\n proc iml;\n mat = {1 2 3, 6 4 5, 7 8 9}; \n tmat = mat`; /* transpose */\n determinant = det(mat); /* get the determinant */\n diagonal_vector = vecdiag(mat); /* get the diagonal as a vector */\n diagonal_mat = diag(mat); /* get the diagonal as a square matrix */\n /* with 0 on off-diagonal entries */\n \n dm = diag({1 2 3}); /* make a square matrix with vector as the diagonal */\n \n print tmat, determinant, diagonal_vector, diagonal_mat, dm;\n quit;\n:::\n:::\n\n\n::: callout-demo\n\n### Matrix Inverse\n\nThe other important matrix-related function is the inverse. \nIn R, `A^-1` will get you the elementwise reciprocal of the matrix. \nNot exactly what we'd like to see... Instead, in R and SAS, we use the `solve()` function. \nThe inverse is defined as the matrix B such that `AB = I` where `I` is the identity matrix (1's on diagonal, 0's off-diagonal). \nSo if we `solve(A)` (in R) or `solve(A, diag(n))` in SAS (where n is a vector of 1s the size of A), we will get the inverse matrix. \nIn Python, we use the `np.linalg.inv()` function to invert a matrix, which may be a bit more linguistically familiar.\n\n::: panel-tabset\n#### R {.unnumbered}\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmat <- matrix(c(1, 2, 3, 6, 4, 5, 7, 8, 9), nrow = 3, byrow = T)\n\nminv <- solve(mat) # get the inverse\n\nminv\n## [,1] [,2] [,3]\n## [1,] -0.2222222 0.3333333 -0.1111111\n## [2,] -1.0555556 -0.6666667 0.7222222\n## [3,] 1.1111111 0.3333333 -0.4444444\nmat %*% minv \n## [,1] [,2] [,3]\n## [1,] 1 0 0\n## [2,] 0 1 0\n## [3,] 0 0 1\n```\n:::\n\n\n#### Python {.unnumbered}\n\n\n::: {.cell}\n\n```{.python .cell-code}\nimport numpy as np\nmat = np.array([[1, 2, 3],[6, 4, 5],[7, 8, 9]], dtype = int, order ='C')\n\nminv = np.linalg.inv(mat)\nminv\n## array([[-0.22222222, 0.33333333, -0.11111111],\n## [-1.05555556, -0.66666667, 0.72222222],\n## [ 1.11111111, 0.33333333, -0.44444444]])\nnp.dot(mat, minv)\n## array([[ 1.00000000e+00, 0.00000000e+00, 1.11022302e-16],\n## [-8.88178420e-16, 1.00000000e+00, -5.55111512e-16],\n## [ 0.00000000e+00, 2.22044605e-16, 1.00000000e+00]])\nnp.round(np.dot(mat, minv), 2)\n## array([[ 1., 0., 0.],\n## [-0., 1., -0.],\n## [ 0., 0., 1.]])\n```\n:::\n\n\n#### SAS {.unnumbered}\n\n[Documentation](https://documentation.sas.com/?docsetId=imlug&docsetTarget=imlug_langref_sect208.htm&docsetVersion=14.2&locale=en)\n\n```\n proc iml;\n mat = {1 2 3, 6 4 5, 7 8 9};\n\n mat_inv = solve(mat, diag({1 1 1})); /* get the inverse */\n mat_inv2 = inv(mat); /* less efficient and less accurate */\n print mat_inv, mat_inv2;\n\n id = mat * mat_inv;\n id2 = mat * mat_inv2;\n print id, id2; \n quit;\n```\n\n:::\n\n:::\n\n## References\n",
0 commit comments