For this reason, the relevant kind of structure with which we will imbue $M$ is a **constraint on the rank of $M$**. In particular, for a positive integer $r \\ll n$, we will assume that $\\mathrm{rank}(M) = r$.

\n",
"\n",
"With $M$ being *of low rank*, it follows that we can represent $M$ by \n",
"$$\n",
"M = WH^*, \\qquad W, H \\in \\mathbb{R}^{n\\times r}\n",
"$$\n",
"which is seen clearly from the SVD-like representation of $M$\n",
"$$\n",
"M = U \\Sigma V^* \\qquad U, V \\in \\mathbb{R}^{n\\times r}, \\Sigma \\in \\mathbb{R}^{r\\times r}\n",
"$$\n",
"where $\\Sigma$ is a diagonal matrix whose elements $\\sigma_1 \\geq \\sigma_2 \\geq \\cdots \\geq \\sigma_r$ are the **non-zero singular values** of $M$; and $U$ and $V$ are unitary matrices (*i.e.*, $UU^* \\equiv U^* U \\equiv VV^* \\equiv V^* V \\equiv 1$).\n",
"\n",
"#### Example: Singular Value Decomposition"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "toCOTn4NQTBy",
"outputId": "e3690dd6-20dd-43a2-fb09-ab5883977775"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"(50, 50)\n",
"50\n",
"(50, 50)\n"
]
}
],
"source": [
"# settings\n",
"n, r = (50, 5)\n",
"W = np.random.randn(n,r)\n",
"H = np.random.randn(n,r)\n",
"M = W @ H.T\n",
"\n",
"# svd\n",
"U, Sigma, Vtr = svd(M)\n",
"\n",
"print(U.shape)\n",
"print(Sigma.size)\n",
"print(Vtr.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "PHJ6TvS6QTBz"
},
"source": [
"Are we really getting 50 singular values for `M`?"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "u7zeK1AnQTBz",
"outputId": "834a0ba0-4e92-4d04-c70b-d9fc6aef725d"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"[64.538 57.021 51.371 43.676 36.003 0. 0. 0. 0. 0. ]\n"
]
}
],
"source": [
"print(Sigma[:2*r].round(3))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "mFgmSUn_QTB0"
},
"source": [
"Yes, but `n-r` of them are equal to $0$ (up to numerical precision). This means we can get a more efficient representation of `M` in SVD format. We'll take only the first `r` columns of `U`; the first `r` rows of `Vtr`; and only the positive singular values of `Sigma`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true,
"id": "EujDdNzRQTB0"
},
"outputs": [],
"source": [
"U_r = U[:, :r]\n",
"Vtr_r = Vtr[:r, :]\n",
"Sigma_r = Sigma[:r]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true,
"id": "KjLVf8Q1QTB0"
},
"outputs": [],
"source": [
"M2 = U_r @ np.diag(Sigma_r) @ Vtr_r"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "5Nvaas1YQTB0",
"outputId": "b85cbf99-d8fb-473d-dc89-fd53750e62a7"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"root mean-squared error = 1.49e-15\n"
]
}
],
"source": [
"print('root mean-squared error = {:5.3g}'.format(rmse(M2, M)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "OfEA_6JXQTB1"
},
"source": [
"**Note:** `scipy.linalg.svd` computes the entire SVD of the matrix — if we know that a matrix is of low rank, then we can use the more efficient `randomized_svd` from `sklearn.decomposition.truncated_svd`, which returns only the top $k$ (or bottom $k$) singular vectors, where $k$ is an argument to the function. \n",
"\n",
"### Sampling pattern\n",
"\n",
"The last problem concerns which entries of $M$ we get to observe. For example, can the entries of $M$ be observed deterministically? Randomly (if so, uniform or non-uniform)? \n",
"In an ideal scenario, we'd like the entries to be observed uniformly at random. Not only does this make the math easier, but it gives us good recovery results. We'll illustrate below what can go wrong if entries are observed non-uniformly. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true,
"id": "o5JwUslDQTB1"
},
"outputs": [],
"source": [
"# dimension and rank\n",
"m = n = 100\n",
"r = 10\n",
"# construct low rank components of M\n",
"W = np.random.randn(m, r)\n",
"H = np.random.randn(n, r)\n",
"# introduce a \"bizarre structure\"\n",
"W[:20,:] = 1\n",
"H[-20:,:] = .75\n",
"# compute M\n",
"M = W @ H.T\n",
"\n",
"# get a non-uniform collection of indices\n",
"i_nu, j_nu = aNonuniformSampling(m, n, .1)\n",
"# Include a uniform sampling of indices for comparison\n",
"i_u, j_u = np.where(np.random.rand(*M.shape) < .1)\n",
"\n",
"# Construct the subsampled matrices\n",
"M_nu = np.zeros(M.shape)\n",
"M_nu[i_nu, j_nu] = M[i_nu, j_nu]\n",
"\n",
"M_u = np.zeros(M.shape)\n",
"M_u[i_u, j_u] = M[i_u, j_u]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 268
},
"id": "KS92yW1ZQTB1",
"outputId": "5433c823-a1e4-4c3c-9383-3b8a5e5a18ee"
},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"