{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"colab":{"provenance":[]},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":3004,"databundleVersionId":861823}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"3e6ef378-6978-48c9-81a0-26cb800460e6","cell_type":"markdown","source":"# Principal Components Analysis (PCA)","metadata":{}},{"id":"02cc24c9-260e-4840-8872-723f0adc6f41","cell_type":"markdown","source":"## 1. Frame the Problem","metadata":{"id":"02cc24c9-260e-4840-8872-723f0adc6f41"}},{"id":"3ae06ffb-3c85-4f69-98aa-fd033c0b1222","cell_type":"code","source":"SCORING = \"accuracy\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:21:01.700278Z","iopub.execute_input":"2026-04-08T13:21:01.700614Z","iopub.status.idle":"2026-04-08T13:21:01.704931Z","shell.execute_reply.started":"2026-04-08T13:21:01.700590Z","shell.execute_reply":"2026-04-08T13:21:01.704005Z"}},"outputs":[],"execution_count":null},{"id":"6c792e7d-2eda-4f9c-86e1-73c9c91984ef","cell_type":"markdown","source":"### Exercise: Load the data\n\nLoad the `train.csv` and `test.csv` csv files into Pandas storing the columns efficiently using the uint8 data type. Write the files to disk using the Parquet format. Compare the file size of the original csv files with the size of the created Parquet files. Which is smaller? Why?","metadata":{}},{"id":"M2nXzqRak_xw","cell_type":"code","source":"%%bash\n\nls /kaggle/input/digit-recognizer","metadata":{"id":"M2nXzqRak_xw","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:21:12.628614Z","iopub.execute_input":"2026-04-08T13:21:12.629002Z","iopub.status.idle":"2026-04-08T13:21:12.641438Z","shell.execute_reply.started":"2026-04-08T13:21:12.628967Z","shell.execute_reply":"2026-04-08T13:21:12.640078Z"}},"outputs":[],"execution_count":null},{"id":"56609008-0b36-4f55-ad2a-966e5d950e37","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"6fd8049c-6646-4783-a24c-c01d8138dc1d","cell_type":"markdown","source":"### Solution","metadata":{}},{"id":"c1fded7d-02a7-48c3-99e4-d19d001c8ccf","cell_type":"code","source":"import pathlib\n\nimport numpy as np\nimport pandas as pd\n\n\n_train_parquet_filepath = pathlib.Path(\"train.parquet\")\nif not _train_parquet_filepath.exists():\n\n    TRAIN_DF = pd.read_csv(\n        \"/kaggle/input/digit-recognizer/train.csv\",\n        dtype=np.uint8,\n    )\n    \n    _ = TRAIN_DF.to_parquet(_train_parquet_filepath)\n\nelse:\n    TRAIN_DF = pd.read_parquet(\n        _train_parquet_filepath,\n    )\n\n_test_features_parquet_filepath = pathlib.Path(\"test-features.parquet\")\nif not _test_features_parquet_filepath.exists():\n    TEST_FEATURES_DF = pd.read_csv(\n        \"/kaggle/input/digit-recognizer/test.csv\",\n        dtype=np.uint8,\n    )\n\n    _ = TEST_FEATURES_DF.to_parquet(_test_features_parquet_filepath)\nelse:\n    TEST_FEATURES_DF = pd.read_parquet(\n        _test_features_parquet_filepath,\n    )\n","metadata":{"id":"c1fded7d-02a7-48c3-99e4-d19d001c8ccf","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:22:07.004410Z","iopub.execute_input":"2026-04-08T13:22:07.004784Z","iopub.status.idle":"2026-04-08T13:22:07.364085Z","shell.execute_reply.started":"2026-04-08T13:22:07.004752Z","shell.execute_reply":"2026-04-08T13:22:07.362413Z"}},"outputs":[],"execution_count":null},{"id":"b5dc764b-dcd3-4984-a3fb-9d52d1c4c215","cell_type":"code","source":"%%bash\n\nls -lh /kaggle/input/digit-recognizer/*.csv","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:22:10.479466Z","iopub.execute_input":"2026-04-08T13:22:10.479825Z","iopub.status.idle":"2026-04-08T13:22:10.496687Z","shell.execute_reply.started":"2026-04-08T13:22:10.479800Z","shell.execute_reply":"2026-04-08T13:22:10.495642Z"}},"outputs":[],"execution_count":null},{"id":"6609ef7b-9f5a-4436-adf1-b264cc5516d9","cell_type":"code","source":"%%bash\n\nls -lh /kaggle/working/*.parquet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:22:15.303901Z","iopub.execute_input":"2026-04-08T13:22:15.304239Z","iopub.status.idle":"2026-04-08T13:22:15.316544Z","shell.execute_reply.started":"2026-04-08T13:22:15.304210Z","shell.execute_reply":"2026-04-08T13:22:15.315602Z"}},"outputs":[],"execution_count":null},{"id":"05f4e494-4b1b-4146-a3c8-e2e2b3c19dfb","cell_type":"markdown","source":"## 2. Sample the training data for faster experimentation\n\nUse stratified sampling to randomly sample 10% of the training data for use in exploring embeddings models and tuning hyperparameters. Later we will refit the best model using the entire training set.","metadata":{}},{"id":"458a8272-df69-491c-a0cb-f0193f7ff725","cell_type":"code","source":"from sklearn import model_selection\n\n\nPRNG = np.random.RandomState(42)\n\n_train_sample_parquet_filepath = pathlib.Path(\"train-sample.parquet\")\nif not _train_sample_parquet_filepath.exists():\n\n    _stratify = TRAIN_DF.loc[:, \"label\"]\n    _train_size = 0.1\n    train_df, _ = model_selection.train_test_split(\n        TRAIN_DF,\n        random_state=PRNG,\n        stratify=_stratify,\n        train_size=_train_size,\n    )\n    \n    _ = train_df.to_parquet(_train_sample_parquet_filepath)\n\nelse:\n    train_df = pd.read_parquet(\n        _train_sample_parquet_filepath,\n    )\n\ntrain_features_df = train_df.drop(\"label\", axis=1)\ntrain_label = train_df.loc[:, \"label\"]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:32:05.702201Z","iopub.execute_input":"2026-04-08T13:32:05.702510Z","iopub.status.idle":"2026-04-08T13:32:05.761988Z","shell.execute_reply.started":"2026-04-08T13:32:05.702475Z","shell.execute_reply":"2026-04-08T13:32:05.760933Z"}},"outputs":[],"execution_count":null},{"id":"c61a0629-33fd-4321-9b01-c552196cea5e","cell_type":"code","source":"%%bash\n\nls -lh /kaggle/working/*.parquet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:32:08.119245Z","iopub.execute_input":"2026-04-08T13:32:08.119515Z","iopub.status.idle":"2026-04-08T13:32:08.129824Z","shell.execute_reply.started":"2026-04-08T13:32:08.119493Z","shell.execute_reply":"2026-04-08T13:32:08.128859Z"}},"outputs":[],"execution_count":null},{"id":"a9b241c6-de2e-4c9a-ac4f-6c5a010e2929","cell_type":"markdown","source":"## 3. Feature Embedding using PCA\n","metadata":{"id":"a9b241c6-de2e-4c9a-ac4f-6c5a010e2929"}},{"id":"b075803a-2549-4b56-8021-ba06f6ddc206","cell_type":"markdown","source":"### PCA is Singular Value Decomposition (SVD)!","metadata":{}},{"id":"b5a9dc1c-35c3-4932-9006-fecae097803b","cell_type":"code","source":"from scipy import linalg\n\n\nlinalg.svd?","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:33:54.990428Z","iopub.execute_input":"2026-04-08T13:33:54.990708Z","iopub.status.idle":"2026-04-08T13:33:54.996648Z","shell.execute_reply.started":"2026-04-08T13:33:54.990685Z","shell.execute_reply":"2026-04-08T13:33:54.995587Z"}},"outputs":[],"execution_count":null},{"id":"5df322a4-c1b5-4a63-b616-5c6d747557f4","cell_type":"code","source":"centered_train_features_df = train_features_df - train_features_df.mean(axis=0)\nU, s, Vh = linalg.svd(centered_train_features_df)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:34:46.797179Z","iopub.execute_input":"2026-04-08T13:34:46.797463Z","iopub.status.idle":"2026-04-08T13:34:48.947602Z","shell.execute_reply.started":"2026-04-08T13:34:46.797441Z","shell.execute_reply":"2026-04-08T13:34:48.947036Z"}},"outputs":[],"execution_count":null},{"id":"4c38cf78-a328-461f-98bd-8274096bcb84","cell_type":"code","source":"U.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:34:49.037925Z","iopub.execute_input":"2026-04-08T13:34:49.038216Z","iopub.status.idle":"2026-04-08T13:34:49.046373Z","shell.execute_reply.started":"2026-04-08T13:34:49.038192Z","shell.execute_reply":"2026-04-08T13:34:49.045897Z"}},"outputs":[],"execution_count":null},{"id":"0f52ebf6-0615-4d48-b0a6-d5746aec1fa3","cell_type":"code","source":"s.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:35:20.703590Z","iopub.execute_input":"2026-04-08T13:35:20.703933Z","iopub.status.idle":"2026-04-08T13:35:20.709398Z","shell.execute_reply.started":"2026-04-08T13:35:20.703911Z","shell.execute_reply":"2026-04-08T13:35:20.708637Z"}},"outputs":[],"execution_count":null},{"id":"af3e979d-396b-430c-9faf-41e22e7b2e07","cell_type":"code","source":"Vh.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:36:01.657243Z","iopub.execute_input":"2026-04-08T13:36:01.657522Z","iopub.status.idle":"2026-04-08T13:36:01.663156Z","shell.execute_reply.started":"2026-04-08T13:36:01.657493Z","shell.execute_reply":"2026-04-08T13:36:01.662260Z"}},"outputs":[],"execution_count":null},{"id":"27933a87-217c-44fd-af2d-5345d08d1e1f","cell_type":"code","source":"# rows of Vh form the individual principle components\nc0 = Vh[0]\nc1 = Vh[1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:36:51.183866Z","iopub.execute_input":"2026-04-08T13:36:51.184185Z","iopub.status.idle":"2026-04-08T13:36:51.188991Z","shell.execute_reply.started":"2026-04-08T13:36:51.184161Z","shell.execute_reply":"2026-04-08T13:36:51.187979Z"}},"outputs":[],"execution_count":null},{"id":"e8c9eb18-4b15-41dc-a62b-e107d06ccb29","cell_type":"markdown","source":"The new, lower-dimensional features are created by projecting the original features into the feature space spanned by the d largest principal components.\n","metadata":{}},{"id":"d76f65b3-0b20-4e3d-b5da-cf70d5329fbf","cell_type":"code","source":"s","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:37:10.491591Z","iopub.execute_input":"2026-04-08T13:37:10.491921Z","iopub.status.idle":"2026-04-08T13:37:10.502190Z","shell.execute_reply.started":"2026-04-08T13:37:10.491894Z","shell.execute_reply":"2026-04-08T13:37:10.501286Z"}},"outputs":[],"execution_count":null},{"id":"ad4e155c-345a-411d-8f6c-bb61aaf9c869","cell_type":"code","source":"d = 2\npca_train_features_df = centered_train_features_df @ Vh[:d].T","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:39:01.565008Z","iopub.execute_input":"2026-04-08T13:39:01.565319Z","iopub.status.idle":"2026-04-08T13:39:01.575887Z","shell.execute_reply.started":"2026-04-08T13:39:01.565296Z","shell.execute_reply":"2026-04-08T13:39:01.575249Z"}},"outputs":[],"execution_count":null},{"id":"d7bd4110-d06f-425d-98e4-d6c35fd9f8ac","cell_type":"code","source":"pca_train_features_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:39:02.299104Z","iopub.execute_input":"2026-04-08T13:39:02.299389Z","iopub.status.idle":"2026-04-08T13:39:02.308119Z","shell.execute_reply.started":"2026-04-08T13:39:02.299363Z","shell.execute_reply":"2026-04-08T13:39:02.307230Z"}},"outputs":[],"execution_count":null},{"id":"43fc80e6-5e4b-42dc-8584-7dfdb397cb27","cell_type":"markdown","source":"### PCA using Scikit Learn","metadata":{"id":"43fc80e6-5e4b-42dc-8584-7dfdb397cb27"}},{"id":"d50c298d-1634-4776-a8a2-b1cd9a5f4894","cell_type":"code","source":"from sklearn import decomposition, pipeline, preprocessing","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:39:51.534676Z","iopub.execute_input":"2026-04-08T13:39:51.535278Z","iopub.status.idle":"2026-04-08T13:39:51.539537Z","shell.execute_reply.started":"2026-04-08T13:39:51.535251Z","shell.execute_reply":"2026-04-08T13:39:51.538471Z"}},"outputs":[],"execution_count":null},{"id":"8cbd8e19-05df-4ff8-9796-cb95697bd099","cell_type":"code","source":"decomposition.PCA?","metadata":{"id":"8cbd8e19-05df-4ff8-9796-cb95697bd099","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:39:56.764097Z","iopub.execute_input":"2026-04-08T13:39:56.764442Z","iopub.status.idle":"2026-04-08T13:39:56.772232Z","shell.execute_reply.started":"2026-04-08T13:39:56.764410Z","shell.execute_reply":"2026-04-08T13:39:56.771211Z"}},"outputs":[],"execution_count":null},{"id":"6f0ce9d7-edec-41a6-82a5-effec0b6a288","cell_type":"code","source":"from sklearn import feature_selection\n\n\n_pca_hyperparameters = {\n    \"n_components\": 2,\n    \"svd_solver\": \"full\",\n    \"whiten\": False,\n}\n\nfeature_embedding = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparameters)\n).set_output(transform=\"pandas\")\n\nembedded_train_features_df = feature_embedding.fit_transform(\n    train_features_df\n)","metadata":{"id":"6f0ce9d7-edec-41a6-82a5-effec0b6a288","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:47:44.764916Z","iopub.execute_input":"2026-04-08T13:47:44.765252Z","iopub.status.idle":"2026-04-08T13:47:45.433974Z","shell.execute_reply.started":"2026-04-08T13:47:44.765221Z","shell.execute_reply":"2026-04-08T13:47:45.433443Z"}},"outputs":[],"execution_count":null},{"id":"9544fd74-c980-4a44-bb3c-73ac31f44bd0","cell_type":"code","source":"embedded_train_features_df.shape","metadata":{"id":"9544fd74-c980-4a44-bb3c-73ac31f44bd0","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:47:47.193438Z","iopub.execute_input":"2026-04-08T13:47:47.193756Z","iopub.status.idle":"2026-04-08T13:47:47.200916Z","shell.execute_reply.started":"2026-04-08T13:47:47.193712Z","shell.execute_reply":"2026-04-08T13:47:47.199502Z"}},"outputs":[],"execution_count":null},{"id":"1802e69b-8053-4ddf-8367-8f8394a1fe63","cell_type":"code","source":"embedded_train_features_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:47:48.152000Z","iopub.execute_input":"2026-04-08T13:47:48.152532Z","iopub.status.idle":"2026-04-08T13:47:48.161213Z","shell.execute_reply.started":"2026-04-08T13:47:48.152502Z","shell.execute_reply":"2026-04-08T13:47:48.160553Z"}},"outputs":[],"execution_count":null},{"id":"4085ca29-c57a-4f14-9db7-50682791d1da","cell_type":"code","source":"embedded_train_features_df.mean(axis=0)","metadata":{"id":"4085ca29-c57a-4f14-9db7-50682791d1da","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:47:55.336634Z","iopub.execute_input":"2026-04-08T13:47:55.338143Z","iopub.status.idle":"2026-04-08T13:47:55.346889Z","shell.execute_reply.started":"2026-04-08T13:47:55.338106Z","shell.execute_reply":"2026-04-08T13:47:55.345280Z"}},"outputs":[],"execution_count":null},{"id":"0880fbb7-46b8-4b8e-be84-70acbf3276bd","cell_type":"code","source":"embedded_train_features_df.std(axis=0)","metadata":{"id":"0880fbb7-46b8-4b8e-be84-70acbf3276bd","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:47:58.535096Z","iopub.execute_input":"2026-04-08T13:47:58.535355Z","iopub.status.idle":"2026-04-08T13:47:58.542324Z","shell.execute_reply.started":"2026-04-08T13:47:58.535335Z","shell.execute_reply":"2026-04-08T13:47:58.541466Z"}},"outputs":[],"execution_count":null},{"id":"04c40f83-2528-4de4-a14f-fcb87258a55a","cell_type":"markdown","source":"## PCA for Compression","metadata":{}},{"id":"f21d42dd-41ca-4899-a0ac-237ad917b73b","cell_type":"code","source":"_pca_hyperparameters = {\n    \"n_components\": 256,\n    \"svd_solver\": \"full\",\n}\n\nfeature_embedding = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparameters)\n).set_output(transform=\"pandas\")\n\ncompressed_train_features_df = feature_embedding.fit_transform(\n    train_features_df\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:49:31.443899Z","iopub.execute_input":"2026-04-08T13:49:31.444202Z","iopub.status.idle":"2026-04-08T13:49:32.056674Z","shell.execute_reply.started":"2026-04-08T13:49:31.444177Z","shell.execute_reply":"2026-04-08T13:49:32.056179Z"}},"outputs":[],"execution_count":null},{"id":"ce4ba588-bd6d-4fa4-853a-e3ac0ca7adc8","cell_type":"markdown","source":"A good compression algorithm should have low reconstruction error. We can measure our reconstruction error of our trained PCA model using our training set.","metadata":{}},{"id":"9c3d91d7-7ffa-4e31-ba08-fc230217ef2b","cell_type":"code","source":"from sklearn import metrics\n\n\n_decompressed_train_features = (\n    feature_embedding.inverse_transform(\n        compressed_train_features_df\n    )\n)\n\n_pixel_reconstruction_errors = metrics.root_mean_squared_error(\n    train_features_df,\n    _decompressed_train_features,\n    multioutput=\"raw_values\",\n)\n\npixel_reconstruction_errors = pd.Series(\n    _pixel_reconstruction_errors,\n    index=train_features_df.columns\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:51:40.549541Z","iopub.execute_input":"2026-04-08T13:51:40.549940Z","iopub.status.idle":"2026-04-08T13:51:40.663819Z","shell.execute_reply.started":"2026-04-08T13:51:40.549914Z","shell.execute_reply":"2026-04-08T13:51:40.663120Z"}},"outputs":[],"execution_count":null},{"id":"dac229b1-856d-4a32-8bd2-504894275057","cell_type":"code","source":"pixel_reconstruction_errors.describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:51:46.327802Z","iopub.execute_input":"2026-04-08T13:51:46.328143Z","iopub.status.idle":"2026-04-08T13:51:46.337554Z","shell.execute_reply.started":"2026-04-08T13:51:46.328114Z","shell.execute_reply":"2026-04-08T13:51:46.336840Z"}},"outputs":[],"execution_count":null},{"id":"39ca5a8b-b595-4a03-99d0-1d6f5467e87a","cell_type":"code","source":"import matplotlib.pyplot as plt\n\n\n_errors = (\n    pixel_reconstruction_errors.to_numpy()\n                               .reshape(28, 28)\n)\n\nplt.imshow(_errors)\nplt.colorbar()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:51:55.936638Z","iopub.execute_input":"2026-04-08T13:51:55.936991Z","iopub.status.idle":"2026-04-08T13:51:56.096134Z","shell.execute_reply.started":"2026-04-08T13:51:55.936966Z","shell.execute_reply":"2026-04-08T13:51:56.094671Z"}},"outputs":[],"execution_count":null},{"id":"c52c8704-aa6e-4736-8943-41c03cdac0be","cell_type":"markdown","source":"## Explained Variance Ratio","metadata":{}},{"id":"94ac78bc-b8a4-47ae-aa2d-a86188963a5e","cell_type":"code","source":"explained_variance_ratio = (\n    feature_embedding.named_steps[\"pca\"]\n                     .explained_variance_ratio_\n)\n\nexplained_variance_ratio","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:54:29.526382Z","iopub.execute_input":"2026-04-08T13:54:29.526652Z","iopub.status.idle":"2026-04-08T13:54:29.534867Z","shell.execute_reply.started":"2026-04-08T13:54:29.526629Z","shell.execute_reply":"2026-04-08T13:54:29.533657Z"}},"outputs":[],"execution_count":null},{"id":"2d974af8-4a4e-4c61-998e-1c5d4c197e73","cell_type":"code","source":"np.cumsum(explained_variance_ratio)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T12:37:26.375501Z","iopub.execute_input":"2026-04-08T12:37:26.375781Z","iopub.status.idle":"2026-04-08T12:37:26.395481Z","shell.execute_reply.started":"2026-04-08T12:37:26.375752Z","shell.execute_reply":"2026-04-08T12:37:26.394793Z"}},"outputs":[],"execution_count":null},{"id":"09c890cd-ffb5-4f15-9498-20858c827b83","cell_type":"markdown","source":"## Choosing the Number of Dimensions","metadata":{"id":"09c890cd-ffb5-4f15-9498-20858c827b83"}},{"id":"ce3a88ee-4111-4d40-946e-e58dbb8e01f4","cell_type":"code","source":"_pca_hyperparameters = {\n    \"n_components\": None,\n    \"svd_solver\": \"full\",\n}\n\nfeature_embedding = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparameters)\n).set_output(transform=\"pandas\")\n\nembedded_train_features_df = feature_embedding.fit_transform(\n    train_features_df\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:59:04.889835Z","iopub.execute_input":"2026-04-08T13:59:04.890125Z","iopub.status.idle":"2026-04-08T13:59:05.557224Z","shell.execute_reply.started":"2026-04-08T13:59:04.890100Z","shell.execute_reply":"2026-04-08T13:59:05.556621Z"}},"outputs":[],"execution_count":null},{"id":"00d91f88-a0c3-4f5a-9d77-9d5026acf76a","cell_type":"code","source":"cummulative_explained_variance_ratio = np.cumsum(\n    feature_embedding.named_steps[\"pca\"]\n                     .explained_variance_ratio_\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:59:06.756516Z","iopub.execute_input":"2026-04-08T13:59:06.757167Z","iopub.status.idle":"2026-04-08T13:59:06.761301Z","shell.execute_reply.started":"2026-04-08T13:59:06.757118Z","shell.execute_reply":"2026-04-08T13:59:06.760312Z"}},"outputs":[],"execution_count":null},{"id":"7d00aa45-802f-4381-99cf-44f1e9acca01","cell_type":"code","source":"len(cummulative_explained_variance_ratio)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:59:08.298453Z","iopub.execute_input":"2026-04-08T13:59:08.298704Z","iopub.status.idle":"2026-04-08T13:59:08.304452Z","shell.execute_reply.started":"2026-04-08T13:59:08.298685Z","shell.execute_reply":"2026-04-08T13:59:08.303302Z"}},"outputs":[],"execution_count":null},{"id":"c381dfc0-10d2-42f6-9625-c0a48b92c650","cell_type":"code","source":"import matplotlib.pyplot as plt\n\n\nplt.plot(\n    np.arange(1, len(cummulative_explained_variance_ratio) + 1),\n    100 * cummulative_explained_variance_ratio\n)\nplt.ylabel(\"Explained Variance (%)\", rotation=\"vertical\")\nplt.xlabel(\"Number of Dimensions, $d$\")\nplt.xscale(\"log\")\nplt.title(\"What is the right choice for $d$?\")\nplt.grid()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T13:59:58.291080Z","iopub.execute_input":"2026-04-08T13:59:58.291361Z","iopub.status.idle":"2026-04-08T13:59:58.661105Z","shell.execute_reply.started":"2026-04-08T13:59:58.291339Z","shell.execute_reply":"2026-04-08T13:59:58.660459Z"}},"outputs":[],"execution_count":null},{"id":"8c4a0198-0434-4d49-9533-cbc72a797e3d","cell_type":"markdown","source":"Choose a threshold value for explained variance, and then choose the minimum value for $d$ necessary to explain required variance.","metadata":{}},{"id":"429d7b15-f503-4450-b541-e76e6c2a22bb","cell_type":"code","source":"explained_variance_threshold = 0.95\nd = np.argmax(cummulative_explained_variance_ratio >= explained_variance_threshold) + 1\nprint(f\"Dimensions required to explain {explained_variance_threshold *100}% of variation in the original data: {d}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:00:46.762635Z","iopub.execute_input":"2026-04-08T14:00:46.763028Z","iopub.status.idle":"2026-04-08T14:00:46.767638Z","shell.execute_reply.started":"2026-04-08T14:00:46.762997Z","shell.execute_reply":"2026-04-08T14:00:46.766987Z"}},"outputs":[],"execution_count":null},{"id":"0cd8e1a7-8110-4359-a6a9-7fdcf83f7ccb","cell_type":"markdown","source":"How accurate is a model trained on the training set reduced to this number of dimensions using PCA?","metadata":{}},{"id":"b8e7a364-b59c-46e3-b7fb-bd8cbbd0b361","cell_type":"code","source":"from sklearn import ensemble, model_selection\n\n\nPRNG = np.random.RandomState(42)\n\nextra_trees_classifier = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(n_components=d, svd_solver=\"full\"),\n    ensemble.ExtraTreesClassifier(n_jobs=1, random_state=PRNG)\n)\n\n_cv_scores = model_selection.cross_val_score(\n    extra_trees_classifier,\n    train_features_df,\n    train_label,\n    cv=3,\n    n_jobs=-1,\n    scoring=SCORING\n)\n\nprint(f\"CV accuracy: {np.mean(_cv_scores * 100):.4f}%\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:05:09.762227Z","iopub.execute_input":"2026-04-08T14:05:09.762545Z","iopub.status.idle":"2026-04-08T14:05:13.810226Z","shell.execute_reply.started":"2026-04-08T14:05:09.762520Z","shell.execute_reply":"2026-04-08T14:05:13.808454Z"}},"outputs":[],"execution_count":null},{"id":"d0df34e5-d939-48e8-a1da-b91ecfefe0ce","cell_type":"markdown","source":"Another option would be to tune the value of $d$ in order to optimize overall model performance on our classification task. ","metadata":{}},{"id":"857ec2a8-2957-4148-9267-ce0a832e0416","cell_type":"code","source":"from scipy import stats\n\n\n_pca_hyperparams = {\n    \"random_state\": PRNG\n}\n_extra_trees_classifier_hyperparams = {\n    \"random_state\": PRNG,\n    \"n_jobs\": 1,\n    \"warm_start\": True,\n}\n\nextra_trees_classifier = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparams),\n    ensemble.ExtraTreesClassifier(**_extra_trees_classifier_hyperparams)\n)\n\n_param_distributions = {\n    \"pca__n_components\": stats.beta(a=9, b=1),\n}\n\nrandomized_search_cv = model_selection.RandomizedSearchCV(\n    extra_trees_classifier,\n    _param_distributions,\n    cv=3,\n    n_iter=10,\n    n_jobs=-1,\n    scoring=SCORING\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:13:20.977482Z","iopub.execute_input":"2026-04-08T14:13:20.977853Z","iopub.status.idle":"2026-04-08T14:13:20.985566Z","shell.execute_reply.started":"2026-04-08T14:13:20.977817Z","shell.execute_reply":"2026-04-08T14:13:20.984822Z"}},"outputs":[],"execution_count":null},{"id":"75b41e8e-2bcf-46e8-ba2b-c7ac9468a736","cell_type":"code","source":"_ = randomized_search_cv.fit(train_features_df, train_label)","metadata":{"id":"75b41e8e-2bcf-46e8-ba2b-c7ac9468a736","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:13:23.398287Z","iopub.execute_input":"2026-04-08T14:13:23.399781Z","iopub.status.idle":"2026-04-08T14:13:44.227164Z","shell.execute_reply.started":"2026-04-08T14:13:23.399683Z","shell.execute_reply":"2026-04-08T14:13:44.226458Z"}},"outputs":[],"execution_count":null},{"id":"4911012f-d518-4153-b00f-fe35bed3740c","cell_type":"code","source":"_cv_accuracy = randomized_search_cv.best_score_\nprint(f\"Validation accuracy: {_cv_accuracy * 100:.4f}%\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:13:44.228648Z","iopub.execute_input":"2026-04-08T14:13:44.229067Z","iopub.status.idle":"2026-04-08T14:13:44.234365Z","shell.execute_reply.started":"2026-04-08T14:13:44.229044Z","shell.execute_reply":"2026-04-08T14:13:44.233550Z"}},"outputs":[],"execution_count":null},{"id":"g5jH-razjScg","cell_type":"code","source":"randomized_search_cv.best_params_","metadata":{"id":"g5jH-razjScg","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:13:44.235600Z","iopub.execute_input":"2026-04-08T14:13:44.235904Z","iopub.status.idle":"2026-04-08T14:13:44.257783Z","shell.execute_reply.started":"2026-04-08T14:13:44.235875Z","shell.execute_reply":"2026-04-08T14:13:44.256362Z"}},"outputs":[],"execution_count":null},{"id":"9db64a60-fd53-4072-9cc5-bfa5f16e95ec","cell_type":"code","source":"optimal_n_components = (\n    randomized_search_cv.best_estimator_\n                        .named_steps[\"pca\"]\n                        .n_components_\n)\n\nprint(optimal_n_components)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:14:32.744819Z","iopub.execute_input":"2026-04-08T14:14:32.745147Z","iopub.status.idle":"2026-04-08T14:14:32.751612Z","shell.execute_reply.started":"2026-04-08T14:14:32.745116Z","shell.execute_reply":"2026-04-08T14:14:32.750358Z"}},"outputs":[],"execution_count":null},{"id":"ndTpbRTipTYU","cell_type":"markdown","source":"### Randomized PCA","metadata":{"id":"ndTpbRTipTYU"}},{"id":"f792c574-ebd3-4d5c-9ed3-5b8db288b17f","cell_type":"code","source":"%%timeit -n 3 -r 1\n\n_pca_hyperparams = {\n    \"n_components\": 3,\n    \"svd_solver\": \"randomized\",\n    \"random_state\": PRNG\n}\n_extra_trees_classifier_hyperparams = {\n    \"random_state\": PRNG,\n    \"n_jobs\": -1, \n}\n\nextra_trees_classifier = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparams),\n    ensemble.ExtraTreesClassifier(**_extra_trees_classifier_hyperparams)\n)\n\n_ = extra_trees_classifier.fit(train_features_df, train_label)","metadata":{"id":"f792c574-ebd3-4d5c-9ed3-5b8db288b17f","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:20:13.704616Z","iopub.execute_input":"2026-04-08T14:20:13.704928Z","iopub.status.idle":"2026-04-08T14:20:15.377040Z","shell.execute_reply.started":"2026-04-08T14:20:13.704904Z","shell.execute_reply":"2026-04-08T14:20:15.375874Z"}},"outputs":[],"execution_count":null},{"id":"42f522b2-cd34-44ab-91d0-a3a2c11056fe","cell_type":"code","source":"%%timeit -n 3 -r 1\n\n_pca_hyperparams = {\n    \"n_components\": 3,\n    \"svd_solver\": \"full\",\n}\n\n_extra_trees_classifier_hyperparams = {\n    \"random_state\": PRNG,\n    \"n_jobs\": -1, \n}\n\nextra_trees_classifier = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(**_pca_hyperparams),\n    ensemble.ExtraTreesClassifier(**_extra_trees_classifier_hyperparams)\n)\n\n_ = extra_trees_classifier.fit(train_features_df, train_label)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:20:27.119637Z","iopub.execute_input":"2026-04-08T14:20:27.119962Z","iopub.status.idle":"2026-04-08T14:20:29.674232Z","shell.execute_reply.started":"2026-04-08T14:20:27.119939Z","shell.execute_reply":"2026-04-08T14:20:29.673515Z"}},"outputs":[],"execution_count":null},{"id":"A1nT1ZwhpxH8","cell_type":"markdown","source":"### Incremental PCA","metadata":{"id":"A1nT1ZwhpxH8"}},{"id":"3dc1259e-5020-4736-8b34-54e958f6a7d3","cell_type":"code","source":"train_mmap = np.memmap(\n    \"train.mmap\",\n    dtype=np.uint8,\n    mode=\"w+\",\n    shape=(42_000, 785)\n)\n\n_chunksize = 1024\ntrain_chunks = pd.read_csv(\n    \"/kaggle/input/digit-recognizer/train.csv\",\n    chunksize=_chunksize,\n    dtype=np.uint8,\n)\n\nfor i, train_chunk in enumerate(train_chunks):\n    train_mmap[_chunksize * i:_chunksize * (i + 1), :] = train_chunk.to_numpy()\n\n# flush the in-memory changes to disk!\ntrain_mmap.flush()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:23:03.466421Z","iopub.execute_input":"2026-04-08T14:23:03.466749Z","iopub.status.idle":"2026-04-08T14:23:05.218229Z","shell.execute_reply.started":"2026-04-08T14:23:03.466704Z","shell.execute_reply":"2026-04-08T14:23:05.216791Z"}},"outputs":[],"execution_count":null},{"id":"c334679f-f62e-4fe4-8515-e673b002fe1e","cell_type":"code","source":"train_mmap = np.memmap(\n    \"train.mmap\",\n    dtype=np.uint8,\n    mode='r',\n    shape=(42_000, 785)\n)\n\n# split off label from features\ntrain_features_mmap = train_mmap[:, 1:]\ntrain_label_mmap = train_mmap[:, 0]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:24:36.453661Z","iopub.execute_input":"2026-04-08T14:24:36.454030Z","iopub.status.idle":"2026-04-08T14:24:36.460574Z","shell.execute_reply.started":"2026-04-08T14:24:36.454007Z","shell.execute_reply":"2026-04-08T14:24:36.459656Z"}},"outputs":[],"execution_count":null},{"id":"ee1264b5-5ed2-4110-a2d8-a03652dce0f2","cell_type":"code","source":"_incremental_pca_hyperparams = {\n    \"n_components\": optimal_n_components,\n    \"batch_size\": 128,\n}\n_extra_trees_classifier_hyperparams = {\n    \"n_jobs\": 1,\n    \"random_state\": PRNG,\n}\n\nextra_trees_classifier = pipeline.make_pipeline(\n    decomposition.IncrementalPCA(**_incremental_pca_hyperparams),\n    ensemble.ExtraTreesClassifier(**_extra_trees_classifier_hyperparams)\n)\n\n_cv_scores = model_selection.cross_val_score(\n    extra_trees_classifier,\n    train_features_mmap,\n    train_label_mmap,\n    cv=3,\n    n_jobs=-1,\n    scoring=SCORING,\n)\n\nprint(f\"CV Accuracy: {np.mean(_cv_scores)}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:27:29.115384Z","iopub.execute_input":"2026-04-08T14:27:29.116175Z","iopub.status.idle":"2026-04-08T14:27:48.381088Z","shell.execute_reply.started":"2026-04-08T14:27:29.116149Z","shell.execute_reply":"2026-04-08T14:27:48.379083Z"}},"outputs":[],"execution_count":null},{"id":"8c4185d0-9c00-42ec-8550-16dfd0c8b9ae","cell_type":"markdown","source":"## 4. Finetuning PCA Feature Embeddings","metadata":{}},{"id":"6b7cb11c-00c4-4ce7-ac13-fef0ff201097","cell_type":"code","source":"from scipy import stats\n\n\nextra_trees_classifier = pipeline.make_pipeline(\n    preprocessing.StandardScaler(),\n    feature_selection.VarianceThreshold(),\n    decomposition.PCA(),\n    ensemble.ExtraTreesClassifier(\n        n_jobs=1,\n        random_state=PRNG,\n    )\n)\n\n_param_distributions = {\n    \"pca__n_components\": stats.beta(a=9, b=1),\n    \"extratreesclassifier__min_samples_split\": stats.beta(a=1, b=9),\n    \"extratreesclassifier__n_estimators\": np.arange(50, 500, step=1)\n}\n\nrandomized_search_cv = model_selection.RandomizedSearchCV(\n    extra_trees_classifier,\n    _param_distributions,\n    cv=3,\n    n_iter=100,\n    n_jobs=-1,\n    scoring=SCORING,\n    verbose=True,\n)\n\n_ = randomized_search_cv.fit(train_features_df, train_label)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:30:28.505876Z","iopub.execute_input":"2026-04-08T14:30:28.506196Z","iopub.status.idle":"2026-04-08T14:33:15.872851Z","shell.execute_reply.started":"2026-04-08T14:30:28.506174Z","shell.execute_reply":"2026-04-08T14:33:15.872221Z"}},"outputs":[],"execution_count":null},{"id":"1ab258d8-5172-480c-a3d7-11da21096b60","cell_type":"code","source":"print(randomized_search_cv.best_score_)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:33:15.873985Z","iopub.execute_input":"2026-04-08T14:33:15.874198Z","iopub.status.idle":"2026-04-08T14:33:15.879620Z","shell.execute_reply.started":"2026-04-08T14:33:15.874178Z","shell.execute_reply":"2026-04-08T14:33:15.878359Z"}},"outputs":[],"execution_count":null},{"id":"caf7a1e3-e38e-4ac2-98cf-621b631baf3c","cell_type":"code","source":"print(randomized_search_cv.best_params_)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:33:15.880781Z","iopub.execute_input":"2026-04-08T14:33:15.881015Z","iopub.status.idle":"2026-04-08T14:33:15.905161Z","shell.execute_reply.started":"2026-04-08T14:33:15.880994Z","shell.execute_reply":"2026-04-08T14:33:15.904106Z"}},"outputs":[],"execution_count":null},{"id":"04f91860-665a-416f-8a62-85a1aecb6fe7","cell_type":"markdown","source":"## 5. Re-fit the tuned model using the entire training set ","metadata":{}},{"id":"7848923d-4b56-402a-b8ba-0a6f234cd164","cell_type":"code","source":"import joblib\n\n\n_train_features_df = TRAIN_DF.drop(\"label\", axis=1)\n_train_label = TRAIN_DF.loc[:, \"label\"]\n\nbest_estimator = randomized_search_cv.best_estimator_\n_ = best_estimator.fit(\n    _train_features_df,\n    _train_label\n)\n    \nbest_estimator_filepath = pathlib.Path(\"best-estimator.pkl\")\n_ = joblib.dump(\n    best_estimator,\n    best_estimator_filepath\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:33:15.907422Z","iopub.execute_input":"2026-04-08T14:33:15.907837Z","iopub.status.idle":"2026-04-08T14:33:52.877139Z","shell.execute_reply.started":"2026-04-08T14:33:15.907797Z","shell.execute_reply":"2026-04-08T14:33:52.875746Z"}},"outputs":[],"execution_count":null},{"id":"6f4f6e1e-eed6-47ad-93f7-11f443bfe53c","cell_type":"markdown","source":"## Submit predictions","metadata":{}},{"id":"G34kX-qb_Njd","cell_type":"code","source":"_test_predictions = best_estimator.predict(TEST_FEATURES_DF)\n\n_ = (\n    pd.read_csv(\n        \"/kaggle/input/digit-recognizer/sample_submission.csv\",\n        index_col=\"ImageId\"\n    ).assign(\n        Label=_test_predictions\n    ).to_csv(\"submission.csv\", index=True)\n)","metadata":{"id":"G34kX-qb_Njd","trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:33:52.878524Z","iopub.execute_input":"2026-04-08T14:33:52.878826Z","iopub.status.idle":"2026-04-08T14:33:54.929476Z","shell.execute_reply.started":"2026-04-08T14:33:52.878803Z","shell.execute_reply":"2026-04-08T14:33:54.928362Z"}},"outputs":[],"execution_count":null},{"id":"167d5ed4-6a6d-4789-90d1-0d38df91fa93","cell_type":"code","source":"%%bash\n\ncat submission.csv | head -n 5","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-04-08T14:33:54.930489Z","iopub.execute_input":"2026-04-08T14:33:54.930739Z","iopub.status.idle":"2026-04-08T14:33:54.942208Z","shell.execute_reply.started":"2026-04-08T14:33:54.930695Z","shell.execute_reply":"2026-04-08T14:33:54.941199Z"}},"outputs":[],"execution_count":null},{"id":"53c750eb-398c-4002-9b33-d37310421d44","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}