{
"cells": [
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"zero (generic function with 18 methods)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Preliminaries: Teach Julia that functions form a vector space\n",
"import Base.+,Base.*,Base.zero\n",
"+(f::Function,g::Function) = x -> f(x)+g(x)\n",
"*(c::Number,f::Function) = x -> c*f(x)\n",
"*(f::Function,c::Number) = x -> c*f(x)\n",
"zero(Function) = x -> 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Examples of Linear Transformations \n",
"\n",
"Some operations are fairly obviously linear. No basis is needed to see this. It is\n",
"efficient theoretically to treat these operations in one fell swoop in a unified way.\n",
"\n",
"For example the derivative of functions $f(x)$ is obviously linear. Derivatives of sums are sums of derivatives: (f+g)'=f'+g'. Derivatives of constant multiples are constant multiples of derivatives (cf)'=cf'.\n",
"Another function transformation example that is obviously linear is the shift by a constant a: $f(x) \\rightarrow f(x+a):$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"T (generic function with 1 method)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function T(f::Function)\n",
" return x->f(x+1)\n",
"end\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" -0.279415 -0.279415"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[T(sin)(5) sin(6)]"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 2.32168 2.32168"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# An example check that shifting is linear\n",
"# we check at x=5 that 2*T(sin)+3*T(cos) = T(2*sin+3*cos), where T denotes shift by one\n",
"[( 2*T(sin) + 3*T(cos) )(5) T( 2*sin + 3*cos )(5)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another example considers the vector space of $m_1 \\times n_1$ matrices $X$. If we take a constant\n",
"$m_2 \\times m_1$ matrix $B$ and a constant $n_2 \\times n_1$ matrix $A$ then the map $X \\rightarrow BXA^T$ is\n",
"obviously a linear map from a vector space of dimension $m_1n_1$ to a vector space of dimension $m_2n_2$.\n",
"(Check: $ B(c_1 X+c_2 Y)A^T= c_1 BXA^T + c_2 BYA^T$.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1: Derivatives of (a cos x + b sin x) \n",
"Consider the 2 dimensional vector space consisting of linear combinations of \"sine\" and \"cosine\". How can we take the derivative of a function in this vector space?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Derivatives: Method 1, symbolically. Matches the paper and pencil method closely."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"using SymPy"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"x,a,b = Sym.([\"x\",\"a\",\"b\"]);"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$- a \\sin{\\left (x \\right )} + b \\cos{\\left (x \\right )}$$"
],
"text/plain": [
"-a*sin(x) + b*cos(x)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff( a*cos(x) + b*sin(x) ,x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Method 2, matrix-vector. Emphasizes the linear nature of derivatives. (Easy to imagine a numerical implementation.)\n",
"\n",
"$\\begin{pmatrix} a' \\\\b' \\end{pmatrix} = \n",
"\\begin{pmatrix} 0 & 1 \\\\-1 & 0 \\end{pmatrix}\n",
"\\begin{pmatrix} a \\\\b \\end{pmatrix}$\n",
"\n",
"Understanding: $\\begin{pmatrix} a' \\\\b' \\end{pmatrix}$ is shorthand for\n",
"$a\\cos x + b\\sin x$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Example: Differentiate $f(x)=5\\cos x + 2 \\sin x$:\n",
"1. Encode f(x) for computation as the vector $\\begin{pmatrix} 5 \\\\ 2 \\end{pmatrix}.$\n",
"2. Apply $\\begin{pmatrix} 0 & 1 \\\\-1 & 0 \\end{pmatrix}$\n",
"to $\\begin{pmatrix} 5 \\\\ 2 \\end{pmatrix}$\n",
"yielding $\\begin{pmatrix} 2 \\\\ -5 \\end{pmatrix}.$\n",
"3. Decode $\\begin{pmatrix} 2 \\\\ -5 \\end{pmatrix}$ as\n",
"$2 \\cos x -5 \\sin x.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Method 3: no shorthand. Combines method 1 and method 2."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\frac{d}{dx} \\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} a \\\\b \\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} 0 & -1 \\\\1 & 0 \\end{pmatrix}\n",
"\\begin{pmatrix} a \\\\b \\end{pmatrix}\n",
"$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Method 2 is purely numerical. The interpretation is imposed by\n",
"the human. Method 3 can be interpreted as method 2 (matrix times vector) with the labels. \n",
"\n",
"If one associates differently, Method 3 can be interpeted as knowing\n",
"the derivative on the basis functions is sufficient for knowing\n",
"this linear transformation everywhere:\n",
"\n",
"$\\frac{d}{dx} \\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} 0 & -1 \\\\1 & 0 \\end{pmatrix}\n",
"$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Observation: \n",
"Method 1 is straightforward but in the end\n",
"bulky. Method 2 shows that the linear transformation defined\n",
"by differentiating can be encoded as a simple matrix times vector,\n",
"which is very efficient on a computer, and also gets to the algebraic heart of the operation. Method 3 organizes the symbolic with the matrices in a way that points to the generalization.\n",
"
\n",
"Most students of calculus learn that differentiation is linear.\n",
"Derivatives of sums are sums of derivatives ,(f+g)'=f'+g'. Derivatives of\n",
"constant multiples are constant multiples of derivatives (cf)'=cf'."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### With code:"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" -1.29521 -1.29521"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f=([sin cos]*[0 1;-1 0]*[5,2])[1]\n",
"x=rand()\n",
"[f(x) 2sin(x)-5cos(x)] ## Check that it gives the right function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## In general\n",
"\n",
"If $v_1,\\ldots,v_n$ is a basis for a vector space $V$, \n",
"and
$w_1,\\ldots,w_m$ is a basis for a vector space $W$,\n",
"and $T$ is some linear transformation,\n",
"we can write\n",
"\n",
"$$ T[v_1,\\ldots,v_n]\n",
"\\begin{pmatrix} c_1 \\\\ \\vdots \\\\ c_n \\end{pmatrix}\n",
"=\n",
"[w_1,\\ldots,w_m] * A* \\begin{pmatrix} c_1 \\\\ \\vdots \\\\ c_n \\end{pmatrix}$$\n",
"for some $m \\times n$ matrix $A$.\n",
"\n",
"One can associate the equation above concentrating\n",
"on\n",
"$$ T[v_1,\\ldots,v_n]\n",
"=\n",
"[w_1,\\ldots,w_m] * A$$\n",
"\n",
"\n",
"to think of\n",
"$T$ as applied to every basis vector of $V$ to get\n",
"some linear combination of the basis vectors of $W$,\n",
"or one can do \"Method 2\" and think of\n",
"$\\begin{pmatrix} c_1 \\\\ \\vdots \\\\ c_n \\end{pmatrix}$\n",
"as the coefficients in the basis for $V$, and \n",
"$A\\begin{pmatrix} c_1 \\\\ \\vdots \\\\ c_n \\end{pmatrix}$\n",
"as the cooeficients of $Tv$ in the basis for $W$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2: Shifting (a cos x + b sin x)\n",
"\n",
"Convince yourself without matrices that\n",
"$Tf$ defined by $ (Tf)(x)=f(x+\\theta)$ is linear for\n",
"any constant $\\theta$.\n",
"\n",
"With matrices we have\n",
"\n",
"$T \\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} \\cos \\theta & -\\sin \\theta \\\\ \\sin \\theta & \\cos \\theta \\end{pmatrix}\n",
"$\n",
"or\n",
"$T \\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} a \\\\b \\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix} \\sin x& \\cos x \\end{pmatrix}\n",
"\\begin{pmatrix} \\cos \\theta & -\\sin \\theta \\\\ \\sin \\theta & \\cos \\theta \\end{pmatrix}\n",
"\\begin{pmatrix} a \\\\b \\end{pmatrix}\n",
"$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which can be done symbolically but gets a little messier looking. The linear algebra is just tidier."
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$a \\sin{\\left (\\theta + x \\right )} + b \\cos{\\left (\\theta + x \\right )}$$"
],
"text/plain": [
"a*sin(theta + x) + b*cos(theta + x)"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Sym(\"x\")\n",
"f = a*sin(x) + b*cos(x)\n",
"Tf = subs(f,x,x+Sym(\"theta\"))"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$a \\left(\\sin{\\left (\\theta \\right )} \\cos{\\left (x \\right )} + \\sin{\\left (x \\right )} \\cos{\\left (\\theta \\right )}\\right) + b \\left(- \\sin{\\left (\\theta \\right )} \\sin{\\left (x \\right )} + \\cos{\\left (\\theta \\right )} \\cos{\\left (x \\right )}\\right)$$"
],
"text/plain": [
"a*(sin(theta)*cos(x) + sin(x)*cos(theta)) + b*(-sin(theta)*sin(x) + cos(theta)\n",
"*cos(x))"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expand_trig(Tf)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course Example 1 is a special case of Example 2 because\n",
"the derivative is the same as shifting by $\\pi/2$ on this very special vector space."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3. Change of basis for polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose one wants to work with Laguerre polynomials.\n",
"Wikipidia can supply the first few of these for us:\n",
" [Laguerre up to degree 6](https://en.wikipedia.org/wiki/Laguerre_polynomials#The_first_few_polynomials)."
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Rational{Int64},2}:\n",
" 1//1 1//1 1//1 1//1 1//1 \n",
" 0//1 -1//1 -2//1 -3//1 -4//1 \n",
" 0//1 0//1 1//2 3//2 3//1 \n",
" 0//1 0//1 0//1 -1//6 -2//3 \n",
" 0//1 0//1 0//1 0//1 1//24"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = Rational.([\n",
" 1 1 2 6 24\n",
" 0 -1 -4 -18 -96\n",
" 0 0 1 9 72\n",
" 0 0 0 -1 -16\n",
" 0 0 0 0 1])./[1 1 2 6 24]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check from the Wikipedia article that\n",
"$[L_0 \\ L_1 \\ L_2 \\ L_3 \\ L_4]=[1 \\ x \\ x^2 \\ x^3 \\ x^4] * A$"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}1&- x + 1&\\frac{x^{2}}{2} - 2 x + 1&- \\frac{x^{3}}{6} + \\frac{3 x^{2}}{2} - 3 x + 1&\\frac{x^{4}}{24} - \\frac{2 x^{3}}{3} + 3 x^{2} - 4 x + 1\\end{bmatrix}"
],
"text/plain": [
"1×5 Array{SymPy.Sym,2}:\n",
" 1 -x + 1 x^2/2 - 2*x + 1 … x^4/24 - 2*x^3/3 + 3*x^2 - 4*x + 1"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1 x x^2 x^3 x^4]*A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Convince yourself that to obtain\n",
"the coefficients of \n",
"$c_0 L_0 + c_1 L_1 + c_2 L_2 + c_3 L_3$\n",
"in the standard basis $1,x,x^2,x^3$\n",
"one must simply compute\n",
"$A * \\begin{pmatrix} c_0 \\\\ c_1 \\\\ c_2 \\\\ c_3 \\end{pmatrix}.$\n",
"\n",
"
\n",
"\n",
"Notationally, we are saying that
\n",
"$[L_0 \\ L_1 \\ L_2 \\ L_3]*\\begin{pmatrix} c_0 \\\\ c_1 \\\\ c_2 \\\\ c_3 \\end{pmatrix}=[1 \\ x \\ x^2 \\ x^3] * A*\\begin{pmatrix} c_0 \\\\ c_1 \\\\ c_2 \\\\ c_3 \\end{pmatrix}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course inv(A) let's us go the other way\n",
"(from monomial coefficients to Laguerre coefficients)"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Int64,2}:\n",
" 1 1 2 6 24\n",
" 0 -1 -4 -18 -96\n",
" 0 0 2 18 144\n",
" 0 0 0 -6 -96\n",
" 0 0 0 0 24"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Int.(inv(A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus for example
\n",
"$x^3 = 6(L_0 - 3 L_1 + 3 L_2 - 1 L_3).$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: the numbers are pascal's triangle times factorials"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What if we want to differentiate quartics written in a Laguerre polynomial basis but we only know how to differentiate in a monomial basis?\n",
"
\n",
"In the standard basis $1,x,x^2,x^3,x^4$ the derivative is this matrix:"
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Int64,2}:\n",
" 0 1 0 0 0\n",
" 0 0 2 0 0\n",
" 0 0 0 3 0\n",
" 0 0 0 0 4\n",
" 0 0 0 0 0"
]
},
"execution_count": 89,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D=[0 1 0 0 0\n",
" 0 0 2 0 0\n",
" 0 0 0 3 0\n",
" 0 0 0 0 4\n",
" 0 0 0 0 0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" We claim\n",
"$\\frac{d}{dx} [1 \\ x \\ x^2 \\ x^3 \\ x^4] = [1 \\ x \\ x^2 \\ x^3 \\ x^4]*D$"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}1&x&x^{2}&x^{3}&x^{4}\\end{bmatrix}"
],
"text/plain": [
"1×5 Array{SymPy.Sym,2}:\n",
" 1 x x^2 x^3 x^4"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1 x x^2 x^3 x^4]"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{bmatrix}0&1&2 x&3 x^{2}&4 x^{3}\\end{bmatrix}"
],
"text/plain": [
"1×5 Array{SymPy.Sym,2}:\n",
" 0 1 2*x 3*x^2 4*x^3"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[1 x x^2 x^3 x^4]*D"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Int64,2}:\n",
" 0 -1 -1 -1 -1\n",
" 0 0 -1 -1 -1\n",
" 0 0 0 -1 -1\n",
" 0 0 0 0 -1\n",
" 0 0 0 0 0"
]
},
"execution_count": 94,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"## Now the derivative in a Laguerre Basis\n",
"Int.(A\\D*A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's interesting. The pattern seems to suggest that\n",
"
\n",
"$\\frac{d}{dx}L_k(x) = -\\sum_{j=0}^{k-1}L_j(x)$\n",
"which is a true identity.\n",
"\n",
"[The wikipedia article](https://en.wikipedia.org/wiki/Laguerre_polynomials) states\n",
"right after the words Sheffer sequence that\n",
"$$\\frac{d}{dx} L_n = -L_{n-1} + \\frac{d}{dx} L_{n-1}$$\n",
"which you should recognize states that the $n$th column\n",
"of the matrix $A^{-1}DA$ is the same as the $n-1$st column,\n",
"with an extra $-1$ in the $(n-1)$st entry."
]
},
{
"cell_type": "code",
"execution_count": 120,
"metadata": {},
"outputs": [],
"source": [
"# Not working yet\n",
"# SymPy.mpmath[:laguerre](4,0,Sym(\"x\")) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Change of Basis\n",
"The above example shows how the derivative matrix can look when changing basis."
]
},
{
"cell_type": "code",
"execution_count": 106,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Int64,2}:\n",
" 0 1 0 0 0\n",
" 0 0 2 0 0\n",
" 0 0 0 3 0\n",
" 0 0 0 0 4\n",
" 0 0 0 0 0"
]
},
"execution_count": 106,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Derivative in standard basis\n",
"D"
]
},
{
"cell_type": "code",
"execution_count": 107,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Int64,2}:\n",
" 0 -1 -1 -1 -1\n",
" 0 0 -1 -1 -1\n",
" 0 0 0 -1 -1\n",
" 0 0 0 0 -1\n",
" 0 0 0 0 0"
]
},
"execution_count": 107,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Derivative in Laguerre basis\n",
"Int.(A\\D*A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Conclusion: Similar matrices represent the same linear transformation in a different basis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3: Kronecker Products"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's check that the [Kronecker Product](https://en.wikipedia.org/wiki/Kronecker_product) is the matrix\n",
"for the transformation $X\\rightarrow BXA^T$"
]
},
{
"cell_type": "code",
"execution_count": 108,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 2 3 2 4\n",
" 8 7 1 4\n",
" 5 9 7 3\n",
" 1 1 2 6"
]
},
"execution_count": 108,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(1:9,4,4)"
]
},
{
"cell_type": "code",
"execution_count": 109,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 2 9 8 6\n",
" 8 2 2 9\n",
" 2 8 6 6\n",
" 4 1 4 9"
]
},
"execution_count": 109,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = rand(1:9,4,4)"
]
},
{
"cell_type": "code",
"execution_count": 110,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"16×16 Array{Int64,2}:\n",
" 4 18 16 12 6 27 24 18 4 18 16 12 8 36 32 24\n",
" 16 4 4 18 24 6 6 27 16 4 4 18 32 8 8 36\n",
" 4 16 12 12 6 24 18 18 4 16 12 12 8 32 24 24\n",
" 8 2 8 18 12 3 12 27 8 2 8 18 16 4 16 36\n",
" 16 72 64 48 14 63 56 42 2 9 8 6 8 36 32 24\n",
" 64 16 16 72 56 14 14 63 8 2 2 9 32 8 8 36\n",
" 16 64 48 48 14 56 42 42 2 8 6 6 8 32 24 24\n",
" 32 8 32 72 28 7 28 63 4 1 4 9 16 4 16 36\n",
" 10 45 40 30 18 81 72 54 14 63 56 42 6 27 24 18\n",
" 40 10 10 45 72 18 18 81 56 14 14 63 24 6 6 27\n",
" 10 40 30 30 18 72 54 54 14 56 42 42 6 24 18 18\n",
" 20 5 20 45 36 9 36 81 28 7 28 63 12 3 12 27\n",
" 2 9 8 6 2 9 8 6 4 18 16 12 12 54 48 36\n",
" 8 2 2 9 8 2 2 9 16 4 4 18 48 12 12 54\n",
" 2 8 6 6 2 8 6 6 4 16 12 12 12 48 36 36\n",
" 4 1 4 9 4 1 4 9 8 2 8 18 24 6 24 54"
]
},
"execution_count": 110,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kron(A,B)"
]
},
{
"cell_type": "code",
"execution_count": 112,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 6 2 7 2\n",
" 5 2 8 4\n",
" 6 1 4 5\n",
" 4 4 5 3"
]
},
"execution_count": 112,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = rand(1:9,4,4)"
]
},
{
"cell_type": "code",
"execution_count": 115,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"16×2 Array{Int64,2}:\n",
" 1108 1108\n",
" 880 880\n",
" 974 974\n",
" 758 758\n",
" 1950 1950\n",
" 1623 1623\n",
" 1714 1714\n",
" 1395 1395\n",
" 2461 2461\n",
" 2110 2110\n",
" 2186 2186\n",
" 1751 1751\n",
" 1067 1067\n",
" 780 780\n",
" 930 930\n",
" 687 687"
]
},
"execution_count": 115,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# The vec operator strings a matrix column wise into one long column\n",
"[ kron(A,B)*vec(X) vec(B*X*A')]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You might check that kron(A,B) is the matrix of the linear transformation $X \\rightarrow BXA^T$\n",
"in the following basis:"
]
},
{
"cell_type": "code",
"execution_count": 119,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 1 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 1 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 1 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 1 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 1 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 1 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 1 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 1 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 1 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 1 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 1 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 1 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 1\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 1\n",
" 0 0 0 0\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 1\n",
" 0 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4×4 Array{Int64,2}:\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" 0 0 0 1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for j=1:4, i=1:4\n",
" V=zeros(Int,4,4)\n",
" V[i,j]=1\n",
" display(V)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.2",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}