{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Cluster Ensembling via SVD\n\nThere's been some some discussion of ensembling clusters, particularly in light of the sensitivity of results to seeds.\n\n- @thedevastor posted a great notebook; [How to Ensemble Clustering algorithms?](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022/discussion/335078#1843648).\n- @aldparis posted a couple of links [How to Ensemble Clustering Algorithms](https://towardsdatascience.com/how-to-ensemble-clustering-algorithms-bf78d7602265) and \n[Cluster_Ensembles](https://pypi.org/project/Cluster_Ensembles/)\n\nI think @thedevastor's notebook and @aldparis' first link are the same algorithm, known as the __Cluster-based Similarity Partitioning Algorithm (CSPA)__ from a paper by\nStehl & Ghosh. It's a very intuitive algorithm, but unfortunately it's computational and storage complexity are both $O(n^2)$. In our case we have ~100,000 observations, so\n$n^2$ is large.\n\nI looked around for some more ideas; interestingly a lot of the literature is ~20 years old, and there aren't many obviously available packages. [Cluster_Ensembles](https://pypi.org/project/Cluster_Ensembles/) implements __CSPA__, as well as __HyperGraph Partitioning (HGPA)__ and __Meta-Clustering MCLA__, also from the Stehl & Ghosh paper. But this package has some issues (i'll try to do a follow up post on it). There is also an R package implementing these 3 algorithms.\n\nI also found a paper from Boulis and Ostendorf that proposes 3 different algorithms, one of which uses Singular Value Decomposition (SVD). There's a part of the math I'm still trying to fully understand, but the basic idea is:\n- $d$: number of observations\n- $c$: number of clusters\n- $s$: number of systems (author's term; essentially the number of clusterings we wish to ensemble)\n- construct matrix $\\mathbf{R}$ of size $d \\times sc$, where each row contains the cluster posteriors of _all_ systems for the given observation\n- decompose $\\mathbf{R}$ via SVD:\n$$ \\mathbf{R} = \\mathbf{U} \\boldsymbol{\\Sigma} \\mathbf{V}^\\intercal$$\nwhere $\\mathbf{U}$ is $d\\times sc$, $\\boldsymbol{\\Sigma}$ is $sc \\times sc$ and $\\mathbf{V}$ is $sc \\times sc$\n- As with low-rank SVD approximation, let $\\mathbf{\\hat{V}}$ be the first $c$ singular vectors of $\\mathbf{V}$.\n- Create final metaspace $\\mathbf{M}=\\mathbf{R}\\mathbf{\\hat{V}}$, size $d \\times c$. From the paper: 'SVD identifies the most correlated clusters and combines them with linear interpolation.' [This is the part I am still trying to fully understand why it works]\n\nNow we run a clustering on the metaspace $M$ to arrive at a consensus clustering (I just tried k-means). Let's see if it works!\n\nLiterature:\n- [Strehl, Ghosh: Cluster Ensembles: A Knowledge Reuse Framework for Combining Multiple Partitions](http://www.strehl.com/download/strehl-aaai02.pdf)\n- [Constantinos Boulis and Mari Ostendorf: Combining Multiple Clustering Systems](https://link.springer.com/content/pdf/10.1007/978-3-540-30116-5_9.pdf)\n","metadata":{}},{"cell_type":"code","source":"import gc\nfrom pathlib import Path\n\nimport numpy as np \nimport pandas as pd\nfrom scipy.linalg import svd\nimport random\n\nimport seaborn as sns \nimport matplotlib.pyplot as plt\n\nfrom sklearn.mixture import GaussianMixture, BayesianGaussianMixture\nfrom sklearn.cluster import KMeans\nfrom sklearn.preprocessing import MaxAbsScaler, RobustScaler, PowerTransformer, QuantileTransformer\nfrom sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score\nfrom sklearn.datasets import make_blobs \nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.pipeline import Pipeline\n\nimport warnings\nwarnings.simplefilter('ignore')\n\npd.set_option('display.max_columns', None)\npd.set_option('display.float_format', '{:.3f}'.format)\n\nINPUT = Path('../input/tabular-playground-series-jul-2022')\nSEED = 420\nN_COMPONENTS=7","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-26T15:11:22.229879Z","iopub.execute_input":"2022-07-26T15:11:22.230556Z","iopub.status.idle":"2022-07-26T15:11:22.241935Z","shell.execute_reply.started":"2022-07-26T15:11:22.230517Z","shell.execute_reply":"2022-07-26T15:11:22.240807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generate some simple clusters","metadata":{}},{"cell_type":"code","source":"# Just generate 3 clusters\nsample_sizes = [500, 500, 500, ] \nmeans = [[0.0, 0.0], [2.0, 2.0], [-2.0, -2.0]]\nstds = [1.0, 1.0, 1.0]\nn_obs = sum(sample_sizes)\nn_clust = len(sample_sizes)\n\n# Generating dataset\nX, y = make_blobs(n_samples = sample_sizes,        \n                  centers = means,               \n                  n_features = 3, \n                  cluster_std = stds, \n                  random_state = SEED) ","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:22.243878Z","iopub.execute_input":"2022-07-26T15:11:22.244802Z","iopub.status.idle":"2022-07-26T15:11:22.253732Z","shell.execute_reply.started":"2022-07-26T15:11:22.244752Z","shell.execute_reply":"2022-07-26T15:11:22.252670Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(8,8)) \nplt.scatter(X[:, 0], X[:, 1], c=y, s=10, alpha=0.5, cmap='Paired')    \nplt.axis('scaled')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:22.255197Z","iopub.execute_input":"2022-07-26T15:11:22.256176Z","iopub.status.idle":"2022-07-26T15:11:22.426950Z","shell.execute_reply.started":"2022-07-26T15:11:22.256135Z","shell.execute_reply":"2022-07-26T15:11:22.426086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Use GaussianMixture to identify clusters\n- not trying to solve optimally so we get a little variance between the different runs\n- save all the labels and posterior probabilities","metadata":{}},{"cell_type":"code","source":"n_iter=5\nposteriors = np.empty((n_obs, n_iter*n_clust))\nlabels = np.empty((n_obs, n_iter))\nfor i in range(n_iter):\n    seed = i*100\n    gmm = GaussianMixture(n_components=n_clust, n_init=1, random_state=seed)\n    gmm.fit(X)\n    i1 = i*n_clust\n    i2 = i1+n_clust\n    posteriors[:,i1:i2] = gmm.predict_proba(X)\n    labels[:,i] = gmm.predict(X)\n\nprint(posteriors[0])\nprint(labels[0])\n","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:22.428064Z","iopub.execute_input":"2022-07-26T15:11:22.429161Z","iopub.status.idle":"2022-07-26T15:11:23.051085Z","shell.execute_reply.started":"2022-07-26T15:11:22.429100Z","shell.execute_reply":"2022-07-26T15:11:23.050190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Decompose the posterior matrix using the SVD\n- verify that the shapes are what we expect","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:18:59.771929Z","iopub.execute_input":"2022-07-26T15:18:59.772383Z","iopub.status.idle":"2022-07-26T15:18:59.777856Z","shell.execute_reply.started":"2022-07-26T15:18:59.772348Z","shell.execute_reply":"2022-07-26T15:18:59.776553Z"}}},{"cell_type":"code","source":"U, s, Vt = svd(posteriors, full_matrices=False)\nprint(U.shape, s.shape, Vt.shape)","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:23.053710Z","iopub.execute_input":"2022-07-26T15:11:23.054050Z","iopub.status.idle":"2022-07-26T15:11:23.066073Z","shell.execute_reply.started":"2022-07-26T15:11:23.054018Z","shell.execute_reply":"2022-07-26T15:11:23.063614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Use n_clust singular vectors to transform the posteriors\n- visualize the resulting meta space\n- we see that it has 'pushed' the different clusters apart","metadata":{}},{"cell_type":"code","source":"meta = posteriors @ Vt[0:n_clust, :].T\n\nfig = plt.figure(figsize=(12,12))\nax = fig.add_subplot(111, projection = '3d')\n\nax.scatter(meta[:,0],meta[:,1],meta[:,2],s=10, alpha=0.5, c=y, cmap='Paired')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:23.076978Z","iopub.execute_input":"2022-07-26T15:11:23.077946Z","iopub.status.idle":"2022-07-26T15:11:23.363800Z","shell.execute_reply.started":"2022-07-26T15:11:23.077888Z","shell.execute_reply":"2022-07-26T15:11:23.362554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Use a clustering algorithm to extract new clusters\n","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\n\nmeta = MinMaxScaler().fit_transform(meta)\nkmeans = KMeans(n_clusters=n_clust, random_state=SEED).fit(meta)\nconsensus = kmeans.predict(meta)\n\nfig = plt.figure(figsize=(12,12))\nax = fig.add_subplot(111, projection = '3d')\n\nax.scatter(meta[:,0],meta[:,1],meta[:,2],s=10, alpha=0.5, c=consensus, cmap='Paired')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:23.365788Z","iopub.execute_input":"2022-07-26T15:11:23.366221Z","iopub.status.idle":"2022-07-26T15:11:23.681481Z","shell.execute_reply.started":"2022-07-26T15:11:23.366179Z","shell.execute_reply":"2022-07-26T15:11:23.680215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### use adjusted rand score to assess ensemble\n\n- we see a slight improvement vs the original scores\n- note that our initial clusters were relatively easy to identify","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import adjusted_rand_score\n\nfor i in range(n_iter):\n    print(adjusted_rand_score(y, labels[:,i]))\n    \nadjusted_rand_score(y, consensus)","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:23.682958Z","iopub.execute_input":"2022-07-26T15:11:23.683338Z","iopub.status.idle":"2022-07-26T15:11:23.701077Z","shell.execute_reply.started":"2022-07-26T15:11:23.683302Z","shell.execute_reply":"2022-07-26T15:11:23.699859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Try the SVD ensemble on the TPS problem","metadata":{}},{"cell_type":"markdown","source":".20036","metadata":{}},{"cell_type":"code","source":"data = pd.read_csv(INPUT / 'data.csv',\n                   index_col='id')\ndata.info()\n\ncols1 = ['f_07','f_08','f_09','f_10','f_11','f_12','f_13']\ncols2 = ['f_22','f_23','f_24','f_25','f_26','f_27','f_28']\ndf = data[cols1 + cols2]\n\ntransformer = Pipeline(\n    steps = [\n        (\"robust\", RobustScaler()),\n        (\"power\", PowerTransformer()),\n    ]\n)\n\nX = transformer.fit_transform(df)\nX = pd.DataFrame(X, columns = df.columns)","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:23.702740Z","iopub.execute_input":"2022-07-26T15:11:23.703069Z","iopub.status.idle":"2022-07-26T15:11:26.347361Z","shell.execute_reply.started":"2022-07-26T15:11:23.703032Z","shell.execute_reply":"2022-07-26T15:11:26.346234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Try 8 BayesianGaussianMixtures\n- use just 1 n_init and different seeds to introduce a little variance","metadata":{}},{"cell_type":"code","source":"n_iter = 16\nn_clust = 7\n\nposteriors = np.empty((X.shape[0], n_iter*n_clust))\n\nprint('fitting seeds: ',end='')\nfor i in range(n_iter):\n    seed = 100 * i\n    print(seed, end='...')\n    model = BayesianGaussianMixture(n_components=n_clust, \n                                    covariance_type='full',\n                                    max_iter=200,\n                                    n_init=1,\n                                    random_state=seed)\n    idx1 = i*n_clust\n    idx2 = idx1+n_clust\n    posteriors[:,idx1:idx2] = model.fit(X).predict_proba(X)","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:11:26.348768Z","iopub.execute_input":"2022-07-26T15:11:26.349087Z","iopub.status.idle":"2022-07-26T15:16:58.852142Z","shell.execute_reply.started":"2022-07-26T15:11:26.349059Z","shell.execute_reply":"2022-07-26T15:16:58.850731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"U, s, Vt = svd(posteriors, full_matrices=False)\nmeta = posteriors @ Vt[0:n_clust, :].T\nkmeans = KMeans(n_clusters=n_clust, random_state=SEED).fit(meta)\nconsensus = kmeans.predict(meta)","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:16:58.854048Z","iopub.execute_input":"2022-07-26T15:16:58.854845Z","iopub.status.idle":"2022-07-26T15:17:01.179811Z","shell.execute_reply.started":"2022-07-26T15:16:58.854798Z","shell.execute_reply":"2022-07-26T15:17:01.178804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = pd.DataFrame(consensus, index=data.index, columns=['Predicted'])\nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2022-07-26T15:17:01.184393Z","iopub.execute_input":"2022-07-26T15:17:01.185057Z","iopub.status.idle":"2022-07-26T15:17:01.390909Z","shell.execute_reply.started":"2022-07-26T15:17:01.185007Z","shell.execute_reply":"2022-07-26T15:17:01.389659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}