{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Active Subspaces\n",
"\n",
"Sometimes, the behavior of an $N$-dimensional model $f$ can be explained best by a *linear reparameterization* of its inputs variables, i.e. we can write $f(\\mathbf{x}) = g(\\mathbf{y}) = g(\\mathbf{M} \\cdot \\mathbf{x})$ where $\\mathbf{M}$ has size $M \\times N$ and $M < N$. When this happens, we say that $f$ admits an $M$-dimensional *active subspace* with basis given by $\\mathbf{M}$'s rows. Those basis vectors are the main directions of variance of the function $f$.\n",
"\n",
"The main directions are the eigenvectors of the matrix\n",
"\n",
"$\\mathbb{E}[\\nabla f^T \\cdot \\nabla f] = \\begin{pmatrix}\n",
"\\mathbb{E}[f_{x_1} \\cdot f_{x_1}] & \\dots & \\mathbb{E}[f_{x_1} \\cdot f_{x_N}] \\\\\n",
"\\dots & \\dots & \\dots \\\\\n",
"\\mathbb{E}[f_{x_N} \\cdot f_{x_1}] & \\dots & \\mathbb{E}[f_{x_N} \\cdot f_{x_N}]\n",
"\\end{pmatrix}$\n",
"\n",
"whereas the eigenvalues reveal the subspace's dimensionality --that is, a large gap between the $M$-th and $(M+1)$-th eigenvalue indicates that an $M$-dimensional active subspace is present.\n",
"\n",
"The necessary expected values are easy to compute from a tensor decomposition: they are just dot products between tensors. We will show a small demonstration of that in this notebook using a 4D function.\n",
"\n",
"Reference: see e.g. [\"Discovering an Active Subspace in a Single-Diode Solar Cell Model\", P. Constantine et al. (2015)](https://arxiv.org/abs/1406.7607)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import tntorch as tn\n",
"import torch\n",
"torch.set_default_dtype(torch.float64)\n",
"\n",
"def f(X):\n",
" return X[:, 0] * X[:, 1] + X[:, 2]\n",
"\n",
"ticks = 64\n",
"P = 100\n",
"N = 4\n",
"\n",
"X = torch.rand((P, N))\n",
"X *= (ticks-1)\n",
"X = torch.round(X)\n",
"y = f(X)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will fit this function `f` using a low-degree expansion in terms of [Legendre polynomials](pce.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter: 0 | loss: 0.999753 | total time: 0.0681\n",
"iter: 500 | loss: 0.976744 | total time: 0.5568\n",
"iter: 1000 | loss: 0.748542 | total time: 1.0160\n",
"iter: 1500 | loss: 0.136286 | total time: 1.4746\n",
"iter: 2000 | loss: 0.008914 | total time: 1.9377\n",
"iter: 2500 | loss: 0.008340 | total time: 2.3975\n",
"iter: 3000 | loss: 0.007649 | total time: 2.8598\n",
"iter: 3500 | loss: 0.006894 | total time: 3.3183\n",
"iter: 4000 | loss: 0.006212 | total time: 3.7835\n",
"iter: 4203 | loss: 0.006041 | total time: 3.9697 <- converged (tol=0.0001)\n"
]
}
],
"source": [
"t = tn.rand([ticks]*N, ranks_tt=2, ranks_tucker=2, requires_grad=True)\n",
"t.set_factors('legendre')\n",
"\n",
"def loss(t):\n",
" return torch.norm(t[X].torch()-y) / torch.norm(y)\n",
"tn.optimize(t, loss)\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX+klEQVR4nO3debgldX3n8feHRRA3QFpEtnbBENwQWsBlEhyVzQWMGwQjKMo4wkQnyZOg4+M6juiMJsG4YUTgUUGjEhBRaQmCZmRpFFklNAQEbDZBFkUG8Dt/VF043tzld7v73HO67/v1POc5Vb/6narvLTj9OfWrOnVSVUiSNJt1Rl2AJGnNYGBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBha0JIcmOT0UdcxkyTfT/LmUdchrTfqAqT5kOQaYHPggYHmY6vqcOBLIylKWsMYGFpIXl5V3xt1EdKayiEpLWhJDk7yw4H5PZJckeSOJJ9KctbgcFCSNyW5PMntSb6bZNuBZZXkrUmuTPKrJJ9MZ4N+/ukDfRcluSfJ45JskuTUJLf06z01yVbT1Pu+JF8cmF/cb3e9fv4xST6fZEWSG5L8zyTr9sue0v89dyS5NclXVu/e1NrOwJB6STYDvga8E3gscAXwvIHl+wLvAv4EWAT8ADhh0mpeBjwHeCbwWmDPqroX+AZwwEC/1wJnVdXNdO/DLwDbAtsA9wD/sJJ/xrHA/cBTgGcDewATgfdB4HRgE2Ar4BMruQ0tUAaGFpJ/7j/pTzzeMmn5PsClVfWNqrofOAq4cWD5W4EPV9Xl/fL/Bew4eJQBHFlVv6qqnwNnAjv27V8G9h/o96d9G1X1y6r6elX9pqruAj4E/PFc/7gkm/d/wzuq6td9GP3twHbvowulJ1TVb6vqh9OsSpqSgaGFZL+q2njg8blJy58AXDcxU92dOa8fWL4t8PcTgQPcBgTYcqDPYMD8BnhkP30msFGSXZMspguSkwCSbJTks0muTXIncDaw8cRQ0hxsC6wPrBio8bPA4/rlf93Xe16SS5O8aY7r1wLnSW/pISvohmoASJLBebow+VBVzfmqqqp6IMlX6YalbgJO7Y8mAP4S+ANg16q6McmOwE/o/nGf7NfARgPzj59U373AZv0R0OQabgTe0v9tLwC+l+Tsqlo+179HC5NHGNJDvgU8I8l+/Unkw/j9f5A/A7wzydPgwRPMr5nD+r8MvA44sJ+e8Ci68xa/SrIp8N4Z1nEh8EdJtknyGLrzLQBU1Qq6cxQfS/LoJOskeXKSP+7rfc3AyfTbgQJ+N4f6tcAZGFpIvpnk7oHHSYMLq+pW4DXAR4FfAjsAy+g+tVNVJwEfAU7sh44uAfZu3XhVnUt3hPAE4NsDi/4OeDhwK3AO8J0Z1rEU+ApwEXABcOqkLm8AHgZcRhcKXwO26Jc9Bzg3yd3AKcDbq+rq1vql+ANK0tSSrEN3DuPAqjpz1PVIo+YRhjQgyZ5JNk6yAd0ltKH71C8teAaG9PueC1xFNzz0crorq+4ZbUnSeHBISpLUxCMMSVKTtfJ7GJtttlktXrx41GVI0hrlggsuuLWqFk23fK0MjMWLF7Ns2bJRlyFJa5Qk18603CEpSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUpO18pveq2rxEd8adQkjdc2RLx11CZLGkEcYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmQwuMJFsnOTPJZUkuTfL2vn3TJEuTXNk/b9K3J8lRSZYnuSjJTgPrOqjvf2WSg4ZVsyRpesM8wrgf+Muq2gHYDTgsyQ7AEcAZVbUdcEY/D7A3sF3/OBT4NHQBA7wX2BXYBXjvRMhIkubP0AKjqlZU1Y/76buAy4EtgX2B4/puxwH79dP7AsdX5xxg4yRbAHsCS6vqtqq6HVgK7DWsuiVJU5uXcxhJFgPPBs4FNq+qFf2iG4HN++ktgesGXnZ93zZd++RtHJpkWZJlt9xyy+r9AyRJww+MJI8Evg68o6ruHFxWVQXU6thOVR1dVUuqasmiRYtWxyolSQOGGhhJ1qcLiy9V1Tf65pv6oSb655v79huArQdevlXfNl27JGkeDfMqqQCfBy6vqo8PLDoFmLjS6SDg5IH2N/RXS+0G3NEPXX0X2CPJJv3J7j36NknSPFpviOt+PvBnwMVJLuzb3gUcCXw1ySHAtcBr+2WnAfsAy4HfAG8EqKrbknwQOL/v94Gqum2IdUuSpjC0wKiqHwKZZvGLpuhfwGHTrOsY4JjVV50kaa78prckqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmQwuMJMckuTnJJQNt70tyQ5IL+8c+A8vemWR5kiuS7DnQvlfftjzJEcOqV5I0s2EeYRwL7DVF+99W1Y794zSAJDsA+wNP61/zqSTrJlkX+CSwN7ADcEDfV5I0z9Yb1oqr6uwkixu77wucWFX3Av+eZDmwS79seVVdDZDkxL7vZau7XknSzJqOMNJ5fZL39PPbJNllttdN4/AkF/VDVpv0bVsC1w30ub5vm659qhoPTbIsybJbbrllJUuTJE2ndUjqU8BzgQP6+bvohorm6tPAk4EdgRXAx1ZiHVOqqqOraklVLVm0aNHqWq0kqdc6JLVrVe2U5CcAVXV7kofNdWNVddPEdJLPAaf2szcAWw903apvY4Z2SdI8aj3CuK8/AV0ASRYBv5vrxpJsMTD7SmDiCqpTgP2TbJDkicB2wHnA+cB2SZ7YB9T+fV9J0jxrPcI4CjgJeFySDwGvBt490wuSnADsDmyW5HrgvcDuSXakC55rgP8CUFWXJvkq3cns+4HDquqBfj2HA98F1gWOqapL5/D3SZJWk6bAqKovJbkAeBEQYL+qunyW1xwwRfPnZ+j/IeBDU7SfBpzWUqckaXhmDIwkmw7M3gycMLisqm4bVmGSpPEy2xHGBXTDRwG2AW7vpzcGfg48cZjFSZLGx4wnvavqiVX1JOB7wMurarOqeizwMuD0+ShQkjQeWq+S2m3iNh4AVfVt4HnDKUmSNI5ar5L6RZJ3A1/s5w8EfjGckiRJ46j1COMAYBHdpbUnAY/joW99S5IWgNbLam8D3j7kWiRJY6wpMJKcSf8t70FV9Z9Xe0WSpLHUeg7jrwamNwReRfeNbEnSAtE6JHXBpKZ/TXLeEOqRJI2p1iGpwW98rwPsDDxmKBVJksZS65DU4De+7wf+HThkWEVJksZPa2D8YVX9drAhyQZDqEeSNKZav4fxf6do+9HqLESSNN5mu1vt4+l+Q/vhSZ5NNyQF8GhgoyHXJkkaI7MNSe0JHEz306gfH2i/C3jXkGqSJI2hGQOjqo4Djkvyqqr6+jzVJEkaQ7MNSb2+qr4ILE7yF5OXV9XHp3iZJGktNNuQ1CP650cOuxBJ0nibbUjqs/3z++enHEnSuGr9pvci4C3A4sHXVNWbhlOWJGnctH5x72TgB3Q/1frA8MqRJI2r1sDYqKr+ZqiVSJLGWus3vU9Nss9QK5EkjbXWwHg7XWjck+TOJHcluXOYhUmSxkvr72E8atiFSJLGW+tVUjtN0XwHcG1V+ct7krQAtJ70/hSwE3BxP/8M4BLgMUn+a1WdPoziJEnjo/Ucxi+AZ1fVzlW1M7AjcDXwEuCjQ6pNkjRGWgPjqVV16cRMVV0GbF9VVw+nLEnSuGkdkro0yaeBE/v51wGX9b+6d99QKpMkjZXWI4yDgeXAO/rH1X3bfcALV39ZkqRx03pZ7T3Ax/rHZHev1ookSWOp9bLa7YAPAzsAG060V9WThlSXJGnMtA5JfQH4NHA/3RDU8cAXh1WUJGn8tAbGw6vqDCBVdW1VvQ946UwvSHJMkpuTXDLQtmmSpUmu7J836duT5Kgky5NcNPhFwSQH9f2vTHLQ3P9ESdLq0BoY9yZZB7gyyeFJXsnsv8J3LLDXpLYjgDOqajvgjH4eYG9gu/5xKN3RDEk2Bd4L7ArsArx3ImQkSfNrLjcf3Aj4c2Bn4M+AGT/tV9XZwG2TmvcFjuunjwP2G2g/vjrnABsn2QLYE1haVbdV1e3AUv5jCEmS5kHrVVLn95N3A29che1tXlUr+ukbgc376S2B6wb6Xd+3Tdf+HyQ5lO7ohG222WYVSpQkTWXGwEhyykzLq+oVK7vhqqoktbKvn2J9RwNHAyxZsmS1rVeS1JntCOO5dJ/wTwDOBbKK27spyRZVtaIfcrq5b78B2Hqg31Z92w3A7pPav7+KNUiSVsJs5zAeD7wLeDrw93Q3G7y1qs6qqrNWYnun8NC5j4Pofit8ov0N/dVSuwF39ENX3wX2SLJJf7J7j75NkjTPZgyMqnqgqr5TVQcBu9HdHuT7SQ6fbcVJTgB+BPxBkuuTHAIcCbwkyZXAi/t5gNPobjeyHPgc8LZ++7cBHwTO7x8f6NskSfNs1pPe/Q0GXwocACwGjgJOmu11VXXANIteNEXfAg6bZj3HAMfMtj1J0nDNdtL7eLrhqNOA91fVJTP1lyStvWY7wng98Gu672H8efLgOe/QHRg8eoi1SZLGyIyBUVWtX+yTJK3lDARJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU1GEhhJrklycZILkyzr2zZNsjTJlf3zJn17khyVZHmSi5LsNIqaJWmhG+URxguraseqWtLPHwGcUVXbAWf08wB7A9v1j0OBT897pZKksRqS2hc4rp8+DthvoP346pwDbJxkixHUJ0kL2qgCo4DTk1yQ5NC+bfOqWtFP3whs3k9vCVw38Nrr+7bfk+TQJMuSLLvllluGVbckLVjrjWi7L6iqG5I8Dlia5GeDC6uqktRcVlhVRwNHAyxZsmROr5UkzW4kRxhVdUP/fDNwErALcNPEUFP/fHPf/QZg64GXb9W3SZLm0bwHRpJHJHnUxDSwB3AJcApwUN/tIODkfvoU4A391VK7AXcMDF1JkubJKIakNgdOSjKx/S9X1XeSnA98NckhwLXAa/v+pwH7AMuB3wBvnP+SJUnzHhhVdTXwrCnafwm8aIr2Ag6bh9IkSTMYp8tqJUljzMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktRkVHer1Vps8RHfGnUJI3XNkS8ddQnSUHiEIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpqsMYGRZK8kVyRZnuSIUdcjSQvNGhEYSdYFPgnsDewAHJBkh9FWJUkLyxoRGMAuwPKqurqq/h9wIrDviGuSpAVlvVEX0GhL4LqB+euBXQc7JDkUOLSfvTvJFTOsbzPg1tVa4eo10vrykVm7uP9m4P4bOutbNTPVt+1ML1xTAmNWVXU0cHRL3yTLqmrJkEtaada3aqxv1Vjfqlmb61tThqRuALYemN+qb5MkzZM1JTDOB7ZL8sQkDwP2B04ZcU2StKCsEUNSVXV/ksOB7wLrAsdU1aWrsMqmoasRsr5VY32rxvpWzVpbX6pqdRYiSVpLrSlDUpKkETMwJElNFkRgJNk0ydIkV/bPm0zT74EkF/aPoZ5Un+1WJ0k2SPKVfvm5SRYPs56VqO/gJLcM7K83z3N9xyS5Ockl0yxPkqP6+i9KstOY1bd7kjsG9t975rm+rZOcmeSyJJcmefsUfUa2DxvrG9k+TLJhkvOS/LSv7/1T9BnZe7ixvrm/h6tqrX8AHwWO6KePAD4yTb+756medYGrgCcBDwN+Cuwwqc/bgM/00/sDX5nH/dVS38HAP4zwv+kfATsBl0yzfB/g20CA3YBzx6y+3YFTR7j/tgB26qcfBfzbFP+NR7YPG+sb2T7s98kj++n1gXOB3Sb1GeV7uKW+Ob+HF8QRBt1tRI7rp48D9htdKUDbrU4Ga/4a8KIkGaP6RqqqzgZum6HLvsDx1TkH2DjJFvNTXVN9I1VVK6rqx/30XcDldHdUGDSyfdhY38j0++Tufnb9/jH5CqKRvYcb65uzhRIYm1fVin76RmDzafptmGRZknOS7DfEeqa61cnkN8ODfarqfuAO4LFDrGnKbfemqg/gVf1QxdeSbD3F8lFq/RtG6bn9kMG3kzxtVEX0QyXPpvsUOmgs9uEM9cEI92GSdZNcCNwMLK2qafffCN7DLfXBHN/Da01gJPlekkumePzeJ+PqjsWmS9ptq/vK/J8Cf5fkycOuew32TWBxVT0TWMpDn6TU5sd0/789C/gE8M+jKCLJI4GvA++oqjtHUcNMZqlvpPuwqh6oqh3p7jyxS5Knz+f2Z9NQ35zfw2tNYFTVi6vq6VM8TgZumjiU7p9vnmYdN/TPVwPfp/tUMwwttzp5sE+S9YDHAL8cUj2TzVpfVf2yqu7tZ/8R2Hmeams11reTqao7J4YMquo0YP0km81nDUnWp/vH+EtV9Y0puox0H85W3zjsw37bvwLOBPaatGiU7+EHTVffyryH15rAmMUpwEH99EHAyZM7JNkkyQb99GbA84HLhlRPy61OBmt+NfAv/dHRfJi1vklj2a+gG2MeJ6cAb+iv9NkNuGNgWHLkkjx+Yjw7yS5078V5+8ek3/bngcur6uPTdBvZPmypb5T7MMmiJBv30w8HXgL8bFK3kb2HW+pbqffwfJ21H+WDbtzwDOBK4HvApn37EuAf++nnARfTXRF0MXDIkGvah+7Kj6uA/9G3fQB4RT+9IfBPwHLgPOBJ87zPZqvvw8Cl/f46E9h+nus7AVgB3Ec3tn4I8Fbgrf3y0P3o1lX9f88lY1bf4QP77xzgefNc3wvohmYvAi7sH/uMyz5srG9k+xB4JvCTvr5LgPf07WPxHm6sb87vYW8NIklqslCGpCRJq8jAkCQ1MTAkSU0MDElSEwNDktTEwJB48Jr+E5NcleSCJKcleeqo61pZ/Z1cnzfqOrR2MTC04PVf/joJ+H5VPbmqdgbeyfT3HFsT7E733SJptTEwJHghcF9VfWaioap+Cvwwyf/u70l2cZLXwYOf3s9KcnKSq5McmeTA/vcHLp64B1mSY5N8pr+h5b8leVnfvmGSL/R9f5LkhX37wUm+keQ76X675aMT9STZI8mPkvw4yT/191giyTVJ3t+3X5xk+/5mfW8F/nu63zn4T/O0H7WWW2/UBUhj4OnABVO0/wmwI/AsYDPg/CRn98ueBfwh3S3Mr6a7Y8Au6X7o578B7+j7Laa7XfyTgTOTPAU4jO4+mM9Isj1w+sDw14509zC7F7giySeAe4B3Ay+uql8n+RvgL+i+tQtwa1XtlORtwF9V1ZuTfIbu913+z6rtGukhBoY0vRcAJ1TVA3Q3sDwLeA5wJ3B+9fdVSnIVcHr/movpjlgmfLWqfgdcmeRqYPt+vZ8AqKqfJbkWmAiMM6rqjn69lwHbAhsDOwD/2t866WHAjwa2MXFjvgvoQk4aCgND6u6n8+o5vubegenfDcz/jt9/X02+985s9+IZXO8D/bpC93sGB8zymon+0lB4DkOCfwE2SHLoREOSZwK/Al6X7odoFtH97Op5c1z3a5Ks05/XeBJwBfAD4MB+O08Ftunbp3MO8Px+OIskj2i4gusuup82lVYbA0MLXnV34Hwl8OL+stpL6e7k+WW6u33+lC5U/rqqbpzj6n9OFzLfprvL6m+BTwHrJLkY+ApwcD30uwRT1XcL3e8vn5DkIrrhqO1n2e43gVd60lurk3erlYYkybHAqVX1tVHXIq0OHmFIkpp4hCFJauIRhiSpiYEhSWpiYEiSmhgYkqQmBoYkqcn/BwpJcmOkGIF0AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"eigvals, eigvecs = tn.active_subspace(t, bounds=None)\n",
"\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.figure()\n",
"plt.bar(range(N), eigvals.detach().numpy())\n",
"plt.title('Eigenvalues')\n",
"plt.xlabel('Component')\n",
"plt.ylabel('Magnitude')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In view of those eigenvalues, we can conclude that the learned model can be written (almost) perfectly in terms of 2 linearly reparameterized variables."
]
}
],
"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.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
