From 1983fc90fe32b4ef3ff2881761ca363a3839ff3b Mon Sep 17 00:00:00 2001 From: root Date: Tue, 16 Jun 2026 04:01:10 +0000 Subject: [PATCH 01/31] Added svd.rs file, where I will code a pure-Rust SVD implementation that does not depend on OpenBLAS --- rust/lance-linalg/src/svd.rs | 1 + 1 file changed, 1 insertion(+) create mode 100644 rust/lance-linalg/src/svd.rs diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs new file mode 100644 index 00000000000..8b137891791 --- /dev/null +++ b/rust/lance-linalg/src/svd.rs @@ -0,0 +1 @@ + From bb307d0646ce92457ea0e7e370a359063d357667 Mon Sep 17 00:00:00 2001 From: root Date: Sat, 20 Jun 2026 03:57:27 +0000 Subject: [PATCH 02/31] Formed the A^T A matrix, and used Jacobi iteration to eigendecompose the A^T A matrix --- benchmarks/full_report/report.ipynb | 126 ++++++---- benchmarks/sift/Results.ipynb | 2 + ci/check_breaking_changes.py | 16 +- ci/setup_version.py | 8 +- memtest/tests/integration_test.rs | 2 +- notebooks/quickstart.ipynb | 31 ++- notebooks/youtube_transcript_search.ipynb | 78 +++--- pre-commit-config.yaml | 11 + rust/lance-linalg/src/svd.rs | 234 ++++++++++++++++++ .../mem_wal/kv/mem_wal_kv_point_lookup.rs | 17 +- .../mem_wal/vector/hnsw/disk_ann_compare.py | 166 ++++++++++--- .../scripts/python_end_to_end.py | 24 +- tall pylance | 173 +++++++++++++ test_data/v0.30.0_pre_created_at/datagen.py | 17 +- test_debug.py | 50 ++-- 15 files changed, 780 insertions(+), 175 deletions(-) create mode 100644 pre-commit-config.yaml create mode 100644 tall pylance diff --git a/benchmarks/full_report/report.ipynb b/benchmarks/full_report/report.ipynb index 11a4924c4c1..70eb810ab26 100644 --- a/benchmarks/full_report/report.ipynb +++ b/benchmarks/full_report/report.ipynb @@ -8,11 +8,12 @@ "outputs": [], "source": [ "import sys\n", + "\n", "sys.dont_write_bytecode = True\n", "\n", "import os\n", "\n", - "module_path = os.path.abspath(os.path.join('.'))\n", + "module_path = os.path.abspath(os.path.join(\".\"))\n", "if module_path not in sys.path:\n", " sys.path.append(module_path)" ] @@ -61,58 +62,91 @@ " gt_in_sample = knn(in_sample, data, metric, 10)\n", "\n", " print(\"generated gt\")\n", - " \n", + "\n", " with tempfile.TemporaryDirectory() as d:\n", " write_lance(d, data)\n", " ds = lance.dataset(d)\n", "\n", - " for q, target in zip(tqdm(in_sample, desc=\"checking brute force\"), gt_in_sample):\n", - " res = ds.to_table(nearest={\n", - " \"column\": \"vec\",\n", - " \"q\": q,\n", - " \"k\": 10,\n", - " \"metric\": metric,\n", - " }, columns=[\"id\"])\n", + " for q, target in zip(\n", + " tqdm(in_sample, desc=\"checking brute force\"), gt_in_sample\n", + " ):\n", + " res = ds.to_table(\n", + " nearest={\n", + " \"column\": \"vec\",\n", + " \"q\": q,\n", + " \"k\": 10,\n", + " \"metric\": metric,\n", + " },\n", + " columns=[\"id\"],\n", + " )\n", " assert len(np.intersect1d(res[\"id\"].to_numpy(), target)) == 10\n", - " \n", - " ds = ds.create_index(\"vec\", \"IVF_PQ\", metric=metric, num_partitions=num_partitions, num_sub_vectors=num_sub_vectors)\n", - " \n", + "\n", + " ds = ds.create_index(\n", + " \"vec\",\n", + " \"IVF_PQ\",\n", + " metric=metric,\n", + " num_partitions=num_partitions,\n", + " num_sub_vectors=num_sub_vectors,\n", + " )\n", + "\n", " recall_data = []\n", " for nprobes in nprobes_list:\n", " for refine_factor in refine_factor_list:\n", " hits = 0\n", " # check that brute force impl is correct\n", - " for q, target in zip(tqdm(query, desc=f\"out of sample, nprobes={nprobes}, refine={refine_factor}\"), gt):\n", - " res = ds.to_table(nearest={\n", - " \"column\": \"vec\",\n", - " \"q\": q,\n", - " \"k\": 10,\n", - " \"nprobes\": nprobes,\n", - " \"refine_factor\": refine_factor,\n", - " }, columns=[\"id\"])[\"id\"].to_numpy()\n", + " for q, target in zip(\n", + " tqdm(\n", + " query,\n", + " desc=f\"out of sample, nprobes={nprobes}, refine={refine_factor}\",\n", + " ),\n", + " gt,\n", + " ):\n", + " res = ds.to_table(\n", + " nearest={\n", + " \"column\": \"vec\",\n", + " \"q\": q,\n", + " \"k\": 10,\n", + " \"nprobes\": nprobes,\n", + " \"refine_factor\": refine_factor,\n", + " },\n", + " columns=[\"id\"],\n", + " )[\"id\"].to_numpy()\n", " hits += len(np.intersect1d(res, target))\n", - " recall_data.append([\n", - " \"out_of_sample\",\n", - " nprobes,\n", - " refine_factor,\n", - " hits / 10 / len(gt),\n", - " ])\n", + " recall_data.append(\n", + " [\n", + " \"out_of_sample\",\n", + " nprobes,\n", + " refine_factor,\n", + " hits / 10 / len(gt),\n", + " ]\n", + " )\n", " # check that brute force impl is correct\n", - " for q, target in zip(tqdm(in_sample, desc=f\"in sample nprobes={nprobes}, refine={refine_factor}\"), gt_in_sample):\n", - " res = ds.to_table(nearest={\n", - " \"column\": \"vec\",\n", - " \"q\": q,\n", - " \"k\": 10,\n", - " \"nprobes\": nprobes,\n", - " \"refine_factor\": refine_factor,\n", - " }, columns=[\"id\"])[\"id\"].to_numpy()\n", + " for q, target in zip(\n", + " tqdm(\n", + " in_sample,\n", + " desc=f\"in sample nprobes={nprobes}, refine={refine_factor}\",\n", + " ),\n", + " gt_in_sample,\n", + " ):\n", + " res = ds.to_table(\n", + " nearest={\n", + " \"column\": \"vec\",\n", + " \"q\": q,\n", + " \"k\": 10,\n", + " \"nprobes\": nprobes,\n", + " \"refine_factor\": refine_factor,\n", + " },\n", + " columns=[\"id\"],\n", + " )[\"id\"].to_numpy()\n", " hits += len(np.intersect1d(res, target))\n", - " recall_data.append([\n", - " \"in_sample\",\n", - " nprobes,\n", - " refine_factor,\n", - " hits / 10 / len(gt_in_sample),\n", - " ])\n", + " recall_data.append(\n", + " [\n", + " \"in_sample\",\n", + " nprobes,\n", + " refine_factor,\n", + " hits / 10 / len(gt_in_sample),\n", + " ]\n", + " )\n", " return recall_data" ] }, @@ -124,15 +158,19 @@ "outputs": [], "source": [ "def make_plot(recall_data):\n", - " df = pd.DataFrame(recall_data, columns=[\"case\", \"nprobes\", \"refine_factor\", \"recall\"])\n", - " \n", + " df = pd.DataFrame(\n", + " recall_data, columns=[\"case\", \"nprobes\", \"refine_factor\", \"recall\"]\n", + " )\n", + "\n", " num_cases = len(df[\"case\"].unique())\n", " (fig, axs) = plt.subplots(1, 2, figsize=(16, 8))\n", - " \n", + "\n", " for case, ax in zip(df[\"case\"].unique(), axs):\n", " current_case = df[df[\"case\"] == case]\n", " sns.heatmap(\n", - " current_case.drop(columns=[\"case\"]).set_index([\"nprobes\", \"refine_factor\"])[\"recall\"].unstack(),\n", + " current_case.drop(columns=[\"case\"])\n", + " .set_index([\"nprobes\", \"refine_factor\"])[\"recall\"]\n", + " .unstack(),\n", " annot=True,\n", " ax=ax,\n", " ).set(title=f\"Recall -- {case}\")" diff --git a/benchmarks/sift/Results.ipynb b/benchmarks/sift/Results.ipynb index 7758587255e..4062e9d45a9 100644 --- a/benchmarks/sift/Results.ipynb +++ b/benchmarks/sift/Results.ipynb @@ -35,6 +35,7 @@ "outputs": [], "source": [ "import pandas as pd\n", + "\n", "df = pd.read_csv(\"query.csv\")" ] }, @@ -46,6 +47,7 @@ "outputs": [], "source": [ "import seaborn as sns\n", + "\n", "sns.set_style(\"darkgrid\")" ] }, diff --git a/ci/check_breaking_changes.py b/ci/check_breaking_changes.py index aa83d1ae7ee..ce569ca637a 100644 --- a/ci/check_breaking_changes.py +++ b/ci/check_breaking_changes.py @@ -4,6 +4,7 @@ Can also be used as a library to detect breaking changes without version validation. """ + import argparse import os import sys @@ -40,10 +41,17 @@ def detect_breaking_changes(repo, base, head): parser = argparse.ArgumentParser() parser.add_argument("base", help="Base commit/tag for comparison") parser.add_argument("head", help="Head commit/tag for comparison") - parser.add_argument("last_stable_version", nargs="?", help="Last stable version (for validation)") - parser.add_argument("current_version", nargs="?", help="Current version (for validation)") - parser.add_argument("--detect-only", action="store_true", - help="Only detect breaking changes, don't validate version") + parser.add_argument( + "last_stable_version", nargs="?", help="Last stable version (for validation)" + ) + parser.add_argument( + "current_version", nargs="?", help="Current version (for validation)" + ) + parser.add_argument( + "--detect-only", + action="store_true", + help="Only detect breaking changes, don't validate version", + ) args = parser.parse_args() repo = Github(os.environ["GITHUB_TOKEN"]).get_repo(os.environ["GITHUB_REPOSITORY"]) diff --git a/ci/setup_version.py b/ci/setup_version.py index 4b248d8e527..b0f6f514d34 100644 --- a/ci/setup_version.py +++ b/ci/setup_version.py @@ -3,7 +3,7 @@ This script is used to set the pre-release version for beta releases (e.g. 0.10.17-beta.1) when the tag indicates a beta release. -With the new automated release process, stable versions are already +With the new automated release process, stable versions are already updated by bump-my-version during the release workflow. """ @@ -38,9 +38,9 @@ def main(): print(f"Setting beta version: {current_version} -> {args.version[1:]}") else: # For stable releases, version should already match - assert ( - parsed_version.release == current_version_parsed.release - ), f"Version mismatch for stable release: {parsed_version.release} != {current_version_parsed.release}" + assert parsed_version.release == current_version_parsed.release, ( + f"Version mismatch for stable release: {parsed_version.release} != {current_version_parsed.release}" + ) with open("python/Cargo.toml", "w") as f: f.writelines(lines) diff --git a/memtest/tests/integration_test.rs b/memtest/tests/integration_test.rs index b83b50cd3d9..a453efd42c3 100644 --- a/memtest/tests/integration_test.rs +++ b/memtest/tests/integration_test.rs @@ -2,7 +2,7 @@ use libc::{c_void, size_t}; use std::ptr; // Import from the library we're testing -use memtest::{memtest_get_stats, memtest_reset_stats, MemtestStats}; +use memtest::{MemtestStats, memtest_get_stats, memtest_reset_stats}; extern "C" { fn malloc(size: size_t) -> *mut c_void; diff --git a/notebooks/quickstart.ipynb b/notebooks/quickstart.ipynb index 5abd79db2d7..2704edfb34d 100644 --- a/notebooks/quickstart.ipynb +++ b/notebooks/quickstart.ipynb @@ -221,7 +221,7 @@ "shutil.rmtree(\"/tmp/test.lance\", ignore_errors=True)\n", "\n", "tbl = pa.Table.from_pandas(df)\n", - "pa.dataset.write_dataset(tbl, \"/tmp/test.parquet\", format='parquet')\n", + "pa.dataset.write_dataset(tbl, \"/tmp/test.parquet\", format=\"parquet\")\n", "\n", "parquet = pa.dataset.dataset(\"/tmp/test.parquet\")\n", "parquet.to_table().to_pandas()" @@ -542,7 +542,7 @@ } ], "source": [ - "lance.dataset('/tmp/test.lance', version=1).to_table().to_pandas()" + "lance.dataset(\"/tmp/test.lance\", version=1).to_table().to_pandas()" ] }, { @@ -600,7 +600,7 @@ } ], "source": [ - "lance.dataset('/tmp/test.lance', version=2).to_table().to_pandas()" + "lance.dataset(\"/tmp/test.lance\", version=2).to_table().to_pandas()" ] }, { @@ -698,7 +698,7 @@ } ], "source": [ - "lance.dataset('/tmp/test.lance', version=\"stable\").to_table().to_pandas()" + "lance.dataset(\"/tmp/test.lance\", version=\"stable\").to_table().to_pandas()" ] }, { @@ -796,11 +796,13 @@ "\n", "with open(\"sift/sift_base.fvecs\", mode=\"rb\") as fobj:\n", " buf = fobj.read()\n", - " data = np.array(struct.unpack(\"<128000000f\", buf[4 : 4 + 4 * 1000000 * 128])).reshape((1000000, 128))\n", + " data = np.array(\n", + " struct.unpack(\"<128000000f\", buf[4 : 4 + 4 * 1000000 * 128])\n", + " ).reshape((1000000, 128))\n", " dd = dict(zip(range(1000000), data))\n", "\n", "table = vec_to_table(dd)\n", - "lance.write_dataset(table, uri, max_rows_per_group=8192, max_rows_per_file=1024*1024)" + "lance.write_dataset(table, uri, max_rows_per_group=8192, max_rows_per_file=1024 * 1024)" ] }, { @@ -864,6 +866,7 @@ ], "source": [ "import duckdb\n", + "\n", "# if this segfaults make sure duckdb v0.7+ is installed\n", "samples = duckdb.query(\"SELECT vector FROM sift1m USING SAMPLE 100\").to_df().vector\n", "samples" @@ -910,10 +913,12 @@ "import time\n", "\n", "start = time.time()\n", - "tbl = sift1m.to_table(columns=[\"id\"], nearest={\"column\": \"vector\", \"q\": samples[0], \"k\": 10})\n", + "tbl = sift1m.to_table(\n", + " columns=[\"id\"], nearest={\"column\": \"vector\", \"q\": samples[0], \"k\": 10}\n", + ")\n", "end = time.time()\n", "\n", - "print(f\"Time(sec): {end-start}\")\n", + "print(f\"Time(sec): {end - start}\")\n", "print(tbl.to_pandas())" ] }, @@ -993,7 +998,7 @@ "\n", "sift1m.create_index(\n", " \"vector\",\n", - " index_type=\"IVF_PQ\", # IVF_PQ, IVF_HNSW_PQ and IVF_HNSW_SQ are supported\n", + " index_type=\"IVF_PQ\", # IVF_PQ, IVF_HNSW_PQ and IVF_HNSW_SQ are supported\n", " num_partitions=256, # IVF\n", " num_sub_vectors=16, # PQ\n", ")" @@ -1068,7 +1073,7 @@ " start = time.time()\n", " tbl = sift1m.to_table(nearest={\"column\": \"vector\", \"q\": q, \"k\": 10})\n", " end = time.time()\n", - " tot += (end - start)\n", + " tot += end - start\n", "\n", "print(f\"Avg(sec): {tot / len(samples)}\")\n", "print(tbl.to_pandas())" @@ -1425,7 +1430,7 @@ "source": [ "tbl = sift1m.to_table()\n", "tbl = tbl.append_column(\"item_id\", pa.array(range(len(tbl))))\n", - "tbl = tbl.append_column(\"revenue\", pa.array((np.random.randn(len(tbl))+5)*1000))\n", + "tbl = tbl.append_column(\"revenue\", pa.array((np.random.randn(len(tbl)) + 5) * 1000))\n", "tbl.to_pandas()" ] }, @@ -1564,7 +1569,9 @@ } ], "source": [ - "sift1m.to_table(columns=[\"revenue\"], nearest={\"column\": \"vector\", \"q\": samples[0], \"k\": 10}).to_pandas()" + "sift1m.to_table(\n", + " columns=[\"revenue\"], nearest={\"column\": \"vector\", \"q\": samples[0], \"k\": 10}\n", + ").to_pandas()" ] } ], diff --git a/notebooks/youtube_transcript_search.ipynb b/notebooks/youtube_transcript_search.ipynb index 42b18d6e557..32f269f02e8 100644 --- a/notebooks/youtube_transcript_search.ipynb +++ b/notebooks/youtube_transcript_search.ipynb @@ -70,7 +70,7 @@ "source": [ "from datasets import load_dataset\n", "\n", - "data = load_dataset('jamescalam/youtube-transcriptions', split='train')\n", + "data = load_dataset(\"jamescalam/youtube-transcriptions\", split=\"train\")\n", "data" ] }, @@ -132,14 +132,20 @@ " text = vid.text.values\n", " time_end = vid[\"end\"].values\n", " contexts = vid.iloc[:-window:stride, :].copy()\n", - " contexts[\"text\"] = [' '.join(text[start_i:start_i+window])\n", - " for start_i in range(0, len(vid)-window, stride)]\n", - " contexts[\"end\"] = [time_end[start_i+window-1]\n", - " for start_i in range(0, len(vid)-window, stride)] \n", + " contexts[\"text\"] = [\n", + " \" \".join(text[start_i : start_i + window])\n", + " for start_i in range(0, len(vid) - window, stride)\n", + " ]\n", + " contexts[\"end\"] = [\n", + " time_end[start_i + window - 1]\n", + " for start_i in range(0, len(vid) - window, stride)\n", + " ]\n", " return contexts\n", + "\n", " # concat result from all videos\n", " return pd.concat([process_video(vid) for _, vid in raw_df.groupby(\"title\")])\n", "\n", + "\n", "df = contextualize(data.to_pandas(), 20, 4)" ] }, @@ -180,7 +186,6 @@ "metadata": {}, "outputs": [], "source": [ - "import functools\n", "import openai\n", "import ratelimiter\n", "from retry import retry\n", @@ -190,12 +195,14 @@ "# API limit at 60/min == 1/sec\n", "limiter = ratelimiter.RateLimiter(max_calls=0.9, period=1.0)\n", "\n", + "\n", "# Get the embedding with retry\n", "@retry(tries=10, delay=1, max_delay=30, backoff=3, jitter=1)\n", - "def embed_func(c): \n", + "def embed_func(c):\n", " rs = openai.Embedding.create(input=c, engine=embed_model)\n", " return [record[\"embedding\"] for record in rs[\"data\"]]\n", "\n", + "\n", "rate_limited = limiter(embed_func)" ] }, @@ -226,15 +233,19 @@ "\n", "openai.api_key = \"sk-...\"\n", "\n", + "\n", "# We request in batches rather than 1 embedding at a time\n", "def to_batches(arr, batch_size):\n", " length = len(arr)\n", + "\n", " def _chunker(arr):\n", " for start_i in range(0, len(df), batch_size):\n", - " yield arr[start_i:start_i+batch_size]\n", + " yield arr[start_i : start_i + batch_size]\n", + "\n", " # add progress meter\n", " yield from tqdm(_chunker(arr), total=math.ceil(length / batch_size))\n", - " \n", + "\n", + "\n", "batch_size = 1000\n", "batches = to_batches(df.text.values.tolist(), batch_size)\n", "embeds = [emb for c in batches for emb in rate_limited(c)]" @@ -280,10 +291,12 @@ } ], "source": [ - "ds = ds.create_index(\"vector\",\n", - " index_type=\"IVF_PQ\", \n", - " num_partitions=64, # IVF\n", - " num_sub_vectors=96) # PQ" + "ds = ds.create_index(\n", + " \"vector\",\n", + " index_type=\"IVF_PQ\",\n", + " num_partitions=64, # IVF\n", + " num_sub_vectors=96,\n", + ") # PQ" ] }, { @@ -304,28 +317,17 @@ "def create_prompt(query, context):\n", " limit = 3750\n", "\n", - " prompt_start = (\n", - " \"Answer the question based on the context below.\\n\\n\"+\n", - " \"Context:\\n\"\n", - " )\n", - " prompt_end = (\n", - " f\"\\n\\nQuestion: {query}\\nAnswer:\"\n", - " )\n", + " prompt_start = \"Answer the question based on the context below.\\n\\n\" + \"Context:\\n\"\n", + " prompt_end = f\"\\n\\nQuestion: {query}\\nAnswer:\"\n", " # append contexts until hitting limit\n", " for i in range(1, len(context)):\n", " if len(\"\\n\\n---\\n\\n\".join(context.text[:i])) >= limit:\n", " prompt = (\n", - " prompt_start +\n", - " \"\\n\\n---\\n\\n\".join(context.text[:i-1]) +\n", - " prompt_end\n", + " prompt_start + \"\\n\\n---\\n\\n\".join(context.text[: i - 1]) + prompt_end\n", " )\n", " break\n", - " elif i == len(context)-1:\n", - " prompt = (\n", - " prompt_start +\n", - " \"\\n\\n---\\n\\n\".join(context.text) +\n", - " prompt_end\n", - " ) \n", + " elif i == len(context) - 1:\n", + " prompt = prompt_start + \"\\n\\n---\\n\\n\".join(context.text) + prompt_end\n", " return prompt" ] }, @@ -350,16 +352,17 @@ "def complete(prompt):\n", " # query text-davinci-003\n", " res = openai.Completion.create(\n", - " engine='text-davinci-003',\n", + " engine=\"text-davinci-003\",\n", " prompt=prompt,\n", " temperature=0,\n", " max_tokens=400,\n", " top_p=1,\n", " frequency_penalty=0,\n", " presence_penalty=0,\n", - " stop=None\n", + " stop=None,\n", " )\n", - " return res['choices'][0]['text'].strip()\n", + " return res[\"choices\"][0][\"text\"].strip()\n", + "\n", "\n", "# check that it works\n", "query = \"who was the 12th person on the moon and when did they land?\"\n", @@ -381,8 +384,9 @@ " \"k\": 3,\n", " \"q\": emb,\n", " \"nprobes\": 20,\n", - " \"refine_factor\": 100\n", - " }).to_pandas()\n", + " \"refine_factor\": 100,\n", + " }\n", + " ).to_pandas()\n", " prompt = create_prompt(question, context)\n", " return complete(prompt), context.reset_index()" ] @@ -443,8 +447,10 @@ } ], "source": [ - "query = (\"Which training method should I use for sentence transformers \"\n", - " \"when I only have pairs of related sentences?\")\n", + "query = (\n", + " \"Which training method should I use for sentence transformers \"\n", + " \"when I only have pairs of related sentences?\"\n", + ")\n", "completion, context = answer(query)\n", "\n", "print(completion)\n", diff --git a/pre-commit-config.yaml b/pre-commit-config.yaml new file mode 100644 index 00000000000..cbc3411268c --- /dev/null +++ b/pre-commit-config.yaml @@ -0,0 +1,11 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v2.3.0 + hooks: + - id: check-yaml + - id: end-of-file-fixer + - id: trailing-whitespace +- repo: https://github.com/psf/black + rev: 22.10.0 + hooks: + - id: black \ No newline at end of file diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 8b137891791..66732ffcd3e 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -1 +1,235 @@ +const EPS: f64 = 2.220_446_049_250_313e-16; +const MAX_JACOBI_SWEEPS: usize = 200; +/// x² — just multiplication, listed for clarity +#[inline] +fn sq(x: f64) -> f64 { x * x } + +/// Raise f64 to an integer power. +fn powi(mut x: f64, mut n: i32) -> f64 { + if n == 0 { return 1.0; } + if n < 0 { x = 1.0 / x; n = -n; } + let mut result = 1.0f64; + let mut base = x; + let mut exp = n as u32; + while exp > 0 { + if exp & 1 == 1 { result *= base; } + base *= base; + exp >>= 1; + } + result +} + +fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { + let mut s = a.to_vec(); + let mut v = eye(n); + + for _ in 0..MAX_JACOBI_SWEEPS { + let mut max_val = 0.0f64; + let mut p = 0; + let mut q = 1; + for i in 0..n { + for j in i + 1..n { + let val = s[i * n + j].abs(); + if val > max_val { + max_val = val; + p = i; + q = j; + } + } + } + if max_val < EPS * 1e4 { + break; + } + + let s_pq = s[p * n + q]; + let diff = s[q * n + q] - s[p * n + p]; + let theta = if diff.abs() < EPS { + std::f64::consts::FRAC_PI_4 + } else { + 0.5 * (2.0 * s_pq / diff).atan() + }; + let (sn, c) = theta.sin_cos(); + + jacobi_rotate(&mut s, n, p, q, c, sn); + apply_givens_right(&mut v, n, p, q, c, sn); + } + + let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); + (eigenvalues, v) +} + +fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + let s_pp = s[p * n + p]; + let s_qq = s[q * n + q]; + let s_pq = s[p * n + q]; + + s[p * n + p] = sq(c) * s_pp - 2.0 * sn * c * s_pq + sq(sn) * s_qq; + s[q * n + q] = sq(sn) * s_pp + 2.0 * sn * c * s_pq + sq(c) * s_qq; + s[p * n + q] = 0.0; + s[q * n + p] = 0.0; + + for r in 0..n { + if r == p || r == q { continue; } + let s_rp = s[r * n + p]; + let s_rq = s[r * n + q]; + let new_rp = c * s_rp - sn * s_rq; + let new_rq = sn * s_rp + c * s_rq; + s[r * n + p] = new_rp; s[p * n + r] = new_rp; + s[r * n + q] = new_rq; s[q * n + r] = new_rq; + } +} + +fn apply_givens_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + for r in 0..n { + let vp = v[r * n + p]; + let vq = v[r * n + q]; + v[r * n + p] = c * vp - sn * vq; + v[r * n + q] = sn * vp + c * vq; + } +} + +// ============================================================ +// Modified Gram-Schmidt orthonormalization +// ============================================================ + +fn gram_schmidt(cols: &mut Vec>) { + let ncols = cols.len(); + for i in 0..ncols { + for j in 0..i { + let dot: f64 = cols[i].iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); + let cj = cols[j].clone(); + for (a, b) in cols[i].iter_mut().zip(cj.iter()) { + *a -= dot * b; + } + } + let norm = (cols[i].iter().map(|&x| sq(x)).sum::()).sqrt(); + if norm > EPS { + for x in cols[i].iter_mut() { *x /= norm; } + } else { + // Replace with an orthonormal basis vector not already spanned + let dim = cols[i].len(); + 'search: for k in 0..dim { + let mut e = vec![0.0f64; dim]; + e[k] = 1.0; + for j in 0..i { + let dot: f64 = e.iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); + let cj = cols[j].clone(); + for (a, b) in e.iter_mut().zip(cj.iter()) { *a -= dot * b; } + } + let n2 = (e.iter().map(|&x| sq(x)).sum::()).sqrt(); + if n2 > EPS { + for x in e.iter_mut() { *x /= n2; } + cols[i] = e; + break 'search; + } + } + } + } +} + +// ============================================================ +// Matrix helpers +// ============================================================ + +/// C = Aᵀ·A where A is m×n → C is n×n +fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { + let mut c = vec![0f64; n * n]; + for i in 0..n { + for l in 0..m { + let a_li = a[l * n + i]; + for j in 0..n { + c[i * n + j] += a_li * a[l * n + j]; + } + } + } + c +} + +/// y = A·x (A is m×n row-major, x length n → y length m) +fn mat_vec_mul(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { + (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() +} + +fn eye(n: usize) -> Vec { + let mut m = vec![0f64; n * n]; + for i in 0..n { m[i * n + i] = 1.0; } + m +} + +pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { + let ata = mat_mul_atb(a, m, n); + let (eigenvalues, v) = jacobi_eigen(&ata, n); + + let k = m.min(n); + let mut sigma: Vec = eigenvalues + .iter() + .map(|&l| (if l > 0.0 { l } else { 0.0 }).sqrt()) + .collect(); + + let mut u_cols: Vec> = Vec::with_capacity(m); + for i in 0..k { + let vi: Vec = (0..n).map(|r| v[r * n + i]).collect(); + let av = mat_vec_mul(a, &vi, m, n); + if sigma[i] > EPS * 10.0 { + u_cols.push(av.iter().map(|&x| x / sigma[i]).collect()); + } else { + sigma[i] = 0.0; + u_cols.push(vec![0.0; m]); + } + } + for _ in k..m { + u_cols.push(vec![0.0; m]); + } + gram_schmidt(&mut u_cols); + + let mut u = vec![0f64; m * m]; + for (ci, col) in u_cols.iter().enumerate() { + for (ri, &value) in col.iter().enumerate() { + u[ri * m + ci] = value; // was: ri * m — m was undefined + } + } + + let mut order: Vec = (0..k).collect(); + for i in 1..k { + let mut j = i; + while j > 0 && sigma[order[j - 1]] < sigma[order[j]] { + order.swap(j - 1, j); + j -= 1; + } + } + + let sigma_sorted: Vec = order.iter().map(|&i| sigma[i]).collect(); // was: .map[( + + let mut u_sorted = vec![0f64; m * m]; + for (new_col, &old_col) in order.iter().enumerate() { + for r in 0..m { + u_sorted[r * m + new_col] = u[r * m + old_col]; // was: r * m (m undefined) + } + } + for col in k..m { + for r in 0..m { + u_sorted[r * m + col] = u[r * m + col]; // was: k..m (m undefined) + } + } + + let mut vt = vec![0f64; n * n]; + for (new_row, &old_col) in order.iter().enumerate() { + for c in 0..n { + vt[new_row * n + c] = v[c * n + old_col]; // was: V_transpose, n undefined + } + } + for row in k..n { + for c in 0..n { + vt[row * n + c] = v[c * n + row]; // was: n undefined + } + } + + (u_sorted, sigma_sorted, vt) +} + +fn main() { + let a = vec![3.0, 2.0, 2.0, 2.0, 3.0, -2.0]; + let (u, s, vt) = svd(&a, 2, 3); + println!("{:?}", s); +} \ No newline at end of file diff --git a/rust/lance/benches/mem_wal/kv/mem_wal_kv_point_lookup.rs b/rust/lance/benches/mem_wal/kv/mem_wal_kv_point_lookup.rs index d99ea97a70a..025a0f27e2c 100644 --- a/rust/lance/benches/mem_wal/kv/mem_wal_kv_point_lookup.rs +++ b/rust/lance/benches/mem_wal/kv/mem_wal_kv_point_lookup.rs @@ -817,8 +817,12 @@ async fn run_lance( }; while lo < part_end { let hi = (lo + args.batch_rows).min(part_end); - let batch = - make_batch(schema.clone(), &insert_order[lo..hi], args.value_size, key_type); + let batch = make_batch( + schema.clone(), + &insert_order[lo..hi], + args.value_size, + key_type, + ); writer.put(vec![batch]).await?; lo = hi; } @@ -1519,7 +1523,14 @@ async fn run_lance_flushed( let peak_rss_mb = sampler.stop(); println!( "[lance] read p50={:.2}us p95={:.2}us p99={:.2}us mean={:.2}us qps_1t={:.0} qps_{}t={:.0} (hits={hits} miss={misses_resolved}) peak_rss={:.0}MB", - stats.p50_us, stats.p95_us, stats.p99_us, stats.mean_us, read_qps_1t, args.threads, read_qps_nt, peak_rss_mb + stats.p50_us, + stats.p95_us, + stats.p99_us, + stats.mean_us, + read_qps_1t, + args.threads, + read_qps_nt, + peak_rss_mb ); Ok(EngineResult { diff --git a/rust/lance/benches/mem_wal/vector/hnsw/disk_ann_compare.py b/rust/lance/benches/mem_wal/vector/hnsw/disk_ann_compare.py index e3b81123ece..1599c2d38ed 100644 --- a/rust/lance/benches/mem_wal/vector/hnsw/disk_ann_compare.py +++ b/rust/lance/benches/mem_wal/vector/hnsw/disk_ann_compare.py @@ -18,7 +18,11 @@ recall@10 vs p50/p99 latency and QPS. The Lance index is served fully cached (large index_cache_size_bytes). """ -import argparse, json, os, time + +import argparse +import json +import os +import time import numpy as np K = 10 @@ -27,7 +31,9 @@ DIM = 1536 EF_SWEEP = [16, 32, 64, 128, 256] HF_TREE = "https://huggingface.co/api/datasets/KShivendu/dbpedia-entities-openai-1M/tree/main/data" -HF_BASE = "https://huggingface.co/datasets/KShivendu/dbpedia-entities-openai-1M/resolve/main/" +HF_BASE = ( + "https://huggingface.co/datasets/KShivendu/dbpedia-entities-openai-1M/resolve/main/" +) def data_dir(base, rows): @@ -42,10 +48,13 @@ def normalize(x): # ---------------- prepare ---------------- def load_corpus(cache_dir, needed): - import requests, pyarrow.parquet as pq + import requests + import pyarrow.parquet as pq + os.makedirs(cache_dir, exist_ok=True) shards = sorted( - e["path"] for e in requests.get(HF_TREE, timeout=60).json() + e["path"] + for e in requests.get(HF_TREE, timeout=60).json() if e["type"] == "file" and e["path"].endswith(".parquet") ) out = np.empty((needed, DIM), dtype=np.float32) @@ -62,7 +71,7 @@ def load_corpus(cache_dir, needed): col = pq.read_table(local, columns=["openai"]).column("openai") arr = np.stack(col.to_pylist()).astype(np.float32) take = min(len(arr), needed - n) - out[n:n + take] = arr[:take] + out[n : n + take] = arr[:take] n += take print(f" shard {os.path.basename(rel)} -> {take} (cum {n})", flush=True) assert n == needed, f"only got {n}/{needed}" @@ -98,10 +107,13 @@ def cmd_prepare(args): rng = np.random.default_rng(SEED) qidx = rng.choice(args.rows, size=NUM_QUERIES, replace=False) queries = corpus[qidx].copy() - print(f"corpus={len(corpus)} queries={len(queries)} dim={DIM}; computing GT...", flush=True) + print( + f"corpus={len(corpus)} queries={len(queries)} dim={DIM}; computing GT...", + flush=True, + ) t = time.perf_counter() gt = numpy_ground_truth(corpus, queries) - print(f" GT in {time.perf_counter()-t:.1f}s", flush=True) + print(f" GT in {time.perf_counter() - t:.1f}s", flush=True) np.save(os.path.join(d, "corpus.npy"), corpus) np.save(os.path.join(d, "queries.npy"), queries) np.save(os.path.join(d, "gt.npy"), gt) @@ -110,7 +122,9 @@ def cmd_prepare(args): # ---------------- shared run helpers ---------------- def recall_at_k(gt, got): - return sum(len(set(g.tolist()) & set(r.tolist())) for g, r in zip(gt, got)) / (len(gt) * K) + return sum(len(set(g.tolist()) & set(r.tolist())) for g, r in zip(gt, got)) / ( + len(gt) * K + ) def latency_qps(query_fn, queries, repeats=3): @@ -133,56 +147,95 @@ def sweep(name, make_q, params, queries, gt): got = np.stack([qf(v) for v in queries]) rec = recall_at_k(gt, got) p50, p99, qps = latency_qps(qf, queries) - rows.append({"param": p, "recall": rec, "p50_us": p50, "p99_us": p99, "qps": qps}) - print(f" {name} param={p} recall={rec:.4f} p50={p50:.0f}us p99={p99:.0f}us qps={qps:.0f}", flush=True) + rows.append( + {"param": p, "recall": rec, "p50_us": p50, "p99_us": p99, "qps": qps} + ) + print( + f" {name} param={p} recall={rec:.4f} p50={p50:.0f}us p99={p99:.0f}us qps={qps:.0f}", + flush=True, + ) return rows # ---------------- systems ---------------- def run_lance(base, rows, corpus, queries, gt): - import lance, pyarrow as pa, shutil + import lance + import pyarrow as pa + import shutil + uri = os.path.join(base, f"lance_{rows}") shutil.rmtree(uri, ignore_errors=True) - vecs = pa.FixedSizeListArray.from_arrays(pa.array(corpus.reshape(-1), type=pa.float32()), DIM) + vecs = pa.FixedSizeListArray.from_arrays( + pa.array(corpus.reshape(-1), type=pa.float32()), DIM + ) tbl = pa.table({"id": pa.array(np.arange(rows, dtype=np.int64)), "vec": vecs}) ds = lance.write_dataset(tbl, uri, mode="overwrite") # The flushed memtable index is a SINGLE-partition HNSW+SQ, so model it with # num_partitions=1 (nprobes=1); ef is the search knob, like DiskANN/FAISS. t = time.perf_counter() - ds.create_index("vec", "IVF_HNSW_SQ", metric="cosine", num_partitions=1, - m=20, ef_construction=150) + ds.create_index( + "vec", + "IVF_HNSW_SQ", + metric="cosine", + num_partitions=1, + m=20, + ef_construction=150, + ) build_s = time.perf_counter() - t ds = lance.dataset(uri, index_cache_size_bytes=48 * 1024**3) def make_q(ef): def q(v): - return ds.to_table(nearest={"column": "vec", "q": v, "k": K, - "nprobes": 1, "ef": ef}, - columns=["id"]).column("id").to_numpy() + return ( + ds.to_table( + nearest={"column": "vec", "q": v, "k": K, "nprobes": 1, "ef": ef}, + columns=["id"], + ) + .column("id") + .to_numpy() + ) + return q - return {"build_s": build_s, "nlist": 1, "sweep": sweep("lance", make_q, None, queries, gt)} + + return { + "build_s": build_s, + "nlist": 1, + "sweep": sweep("lance", make_q, None, queries, gt), + } def run_lance_flushed(base, rows, corpus, queries, gt, lance_path, id_offset, column): # Open a flushed MemTable generation directly from its dataset path and # benchmark its on-disk IVF_HNSW_SQ index (single partition), fully cached. import lance + ds = lance.dataset(lance_path, index_cache_size_bytes=48 * 1024**3) def make_q(ef): def q(v): - ids = ds.to_table(nearest={"column": column, "q": v, "k": K, - "nprobes": 1, "ef": ef}, - columns=["id"]).column("id").to_numpy() + ids = ( + ds.to_table( + nearest={"column": column, "q": v, "k": K, "nprobes": 1, "ef": ef}, + columns=["id"], + ) + .column("id") + .to_numpy() + ) return ids - id_offset # map flushed-gen id -> corpus index + return q - return {"lance_path": lance_path, "id_offset": id_offset, - "sweep": sweep("lance_flushed", make_q, None, queries, gt)} + + return { + "lance_path": lance_path, + "id_offset": id_offset, + "sweep": sweep("lance_flushed", make_q, None, queries, gt), + } def run_faiss(base, rows, corpus, queries, gt): # Full-precision HNSW reference (shows what no quantization buys). import faiss + index = faiss.IndexHNSWFlat(DIM, 32, faiss.METRIC_INNER_PRODUCT) index.hnsw.efConstruction = 200 t = time.perf_counter() @@ -194,16 +247,20 @@ def make_q(ef): def q(v): index.hnsw.efSearch = ef return index.search(v.reshape(1, -1), K)[1][0] + return q + return {"build_s": build_s, "sweep": sweep("faiss", make_q, None, queries, gt)} def run_faiss_sq(base, rows, corpus, queries, gt): # HNSW + 8-bit scalar quantization — apples-to-apples with Lance IVF_HNSW_SQ. import faiss + try: - index = faiss.IndexHNSWSQ(DIM, faiss.ScalarQuantizer.QT_8bit, 32, - faiss.METRIC_INNER_PRODUCT) + index = faiss.IndexHNSWSQ( + DIM, faiss.ScalarQuantizer.QT_8bit, 32, faiss.METRIC_INNER_PRODUCT + ) except Exception: # Fall back to L2; on unit-normalized vectors L2 ranking == cosine. index = faiss.IndexHNSWSQ(DIM, faiss.ScalarQuantizer.QT_8bit, 32) @@ -218,28 +275,44 @@ def make_q(ef): def q(v): index.hnsw.efSearch = ef return index.search(v.reshape(1, -1), K)[1][0] + return q + return {"build_s": build_s, "sweep": sweep("faiss_sq", make_q, None, queries, gt)} def run_diskann(base, rows, corpus, queries, gt): import diskannpy as dap + idx_dir = os.path.join(base, f"diskann_{rows}") os.makedirs(idx_dir, exist_ok=True) t = time.perf_counter() dap.build_memory_index( - data=corpus, distance_metric="cosine", index_directory=idx_dir, - index_prefix="ann", complexity=150, graph_degree=64, - num_threads=0, alpha=1.2, use_pq_build=False, num_pq_bytes=0, + data=corpus, + distance_metric="cosine", + index_directory=idx_dir, + index_prefix="ann", + complexity=150, + graph_degree=64, + num_threads=0, + alpha=1.2, + use_pq_build=False, + num_pq_bytes=0, ) build_s = time.perf_counter() - t - idx = dap.StaticMemoryIndex(index_directory=idx_dir, index_prefix="ann", - num_threads=0, initial_search_complexity=256) + idx = dap.StaticMemoryIndex( + index_directory=idx_dir, + index_prefix="ann", + num_threads=0, + initial_search_complexity=256, + ) def make_q(L): def q(v): return idx.search(v, k_neighbors=K, complexity=max(L, K)).identifiers + return q + return {"build_s": build_s, "sweep": sweep("diskann", make_q, None, queries, gt)} @@ -250,11 +323,23 @@ def cmd_run(args): gt = np.load(os.path.join(d, "gt.npy")) print(f"=== {args.system} rows={args.rows} corpus={len(corpus)} ===", flush=True) if args.system == "lance_flushed": - res = run_lance_flushed(args.base, args.rows, corpus, queries, gt, - args.lance_path, args.id_offset, args.column) + res = run_lance_flushed( + args.base, + args.rows, + corpus, + queries, + gt, + args.lance_path, + args.id_offset, + args.column, + ) else: - fn = {"lance": run_lance, "faiss": run_faiss, "faiss_sq": run_faiss_sq, - "diskann": run_diskann}[args.system] + fn = { + "lance": run_lance, + "faiss": run_faiss, + "faiss_sq": run_faiss_sq, + "diskann": run_diskann, + }[args.system] res = fn(args.base, args.rows, corpus, queries, gt) res["rows"] = args.rows res["system"] = args.system @@ -267,9 +352,16 @@ def cmd_run(args): def main(): ap = argparse.ArgumentParser() sub = ap.add_subparsers(dest="cmd", required=True) - p = sub.add_parser("prepare"); p.add_argument("--rows", type=int, required=True); p.add_argument("--base", required=True) - r = sub.add_parser("run"); r.add_argument("--rows", type=int, required=True); r.add_argument("--base", required=True); r.add_argument("--system", required=True) - r.add_argument("--lance-path", default=None); r.add_argument("--id-offset", type=int, default=0); r.add_argument("--column", default="vector") + p = sub.add_parser("prepare") + p.add_argument("--rows", type=int, required=True) + p.add_argument("--base", required=True) + r = sub.add_parser("run") + r.add_argument("--rows", type=int, required=True) + r.add_argument("--base", required=True) + r.add_argument("--system", required=True) + r.add_argument("--lance-path", default=None) + r.add_argument("--id-offset", type=int, default=0) + r.add_argument("--column", default="vector") args = ap.parse_args() (cmd_prepare if args.cmd == "prepare" else cmd_run)(args) diff --git a/skills/lance-user-guide/scripts/python_end_to_end.py b/skills/lance-user-guide/scripts/python_end_to_end.py index ec2d02713c9..b366fe93151 100644 --- a/skills/lance-user-guide/scripts/python_end_to_end.py +++ b/skills/lance-user-guide/scripts/python_end_to_end.py @@ -11,16 +11,24 @@ import lance -def _build_fixed_size_vectors(num_rows: int, dim: int) -> tuple[pa.FixedSizeListArray, np.ndarray]: +def _build_fixed_size_vectors( + num_rows: int, dim: int +) -> tuple[pa.FixedSizeListArray, np.ndarray]: vectors = np.random.rand(num_rows, dim).astype("float32") flat = pa.array(vectors.reshape(-1), type=pa.float32()) return pa.FixedSizeListArray.from_arrays(flat, dim), vectors def main() -> None: - parser = argparse.ArgumentParser(description="Minimal Lance write/index/query example") - parser.add_argument("--uri", default="example.lance", help="Dataset URI (directory)") - parser.add_argument("--mode", default="overwrite", choices=["create", "append", "overwrite"]) + parser = argparse.ArgumentParser( + description="Minimal Lance write/index/query example" + ) + parser.add_argument( + "--uri", default="example.lance", help="Dataset URI (directory)" + ) + parser.add_argument( + "--mode", default="overwrite", choices=["create", "append", "overwrite"] + ) parser.add_argument("--rows", type=int, default=1000) parser.add_argument("--dim", type=int, default=32) @@ -40,7 +48,13 @@ def main() -> None: uri = str(Path(args.uri)) vec_arr, vec_np = _build_fixed_size_vectors(args.rows, args.dim) categories = pa.array(["a" if i % 2 == 0 else "b" for i in range(args.rows)]) - table = pa.table({"id": pa.array(range(args.rows), pa.int64()), "category": categories, "vector": vec_arr}) + table = pa.table( + { + "id": pa.array(range(args.rows), pa.int64()), + "category": categories, + "vector": vec_arr, + } + ) ds = lance.write_dataset(table, uri, mode=args.mode) ds = lance.dataset(uri) diff --git a/tall pylance b/tall pylance new file mode 100644 index 00000000000..933b1652748 --- /dev/null +++ b/tall pylance @@ -0,0 +1,173 @@ +diff --git a/benchmarks/full_report/report.ipynb b/benchmarks/full_report/report.ipynb +index 11a4924c..70eb810a 100644 +--- a/benchmarks/full_report/report.ipynb ++++ b/benchmarks/full_report/report.ipynb +@@ -8,11 +8,12 @@ + "outputs": [], + "source": [ + "import sys\n", ++ "\n", + "sys.dont_write_bytecode = True\n", + "\n", + "import os\n", + "\n", +- "module_path = os.path.abspath(os.path.join('.'))\n", ++ "module_path = os.path.abspath(os.path.join(\".\"))\n", + "if module_path not in sys.path:\n", + " sys.path.append(module_path)" + ] +@@ -61,58 +62,91 @@ + " gt_in_sample = knn(in_sample, data, metric, 10)\n", + "\n", + " print(\"generated gt\")\n", +- " \n", ++ "\n", + " with tempfile.TemporaryDirectory() as d:\n", + " write_lance(d, data)\n", + " ds = lance.dataset(d)\n", + "\n", +- " for q, target in zip(tqdm(in_sample, desc=\"checking brute force\"), gt_in_sample):\n", +- " res = ds.to_table(nearest={\n", +- " \"column\": \"vec\",\n", +- " \"q\": q,\n", +- " \"k\": 10,\n", +- " \"metric\": metric,\n", +- " }, columns=[\"id\"])\n", ++ " for q, target in zip(\n", ++ " tqdm(in_sample, desc=\"checking brute force\"), gt_in_sample\n", ++ " ):\n", ++ " res = ds.to_table(\n", ++ " nearest={\n", ++ " \"column\": \"vec\",\n", ++ " \"q\": q,\n", ++ " \"k\": 10,\n", ++ " \"metric\": metric,\n", ++ " },\n", ++ " columns=[\"id\"],\n", ++ " )\n", + " assert len(np.intersect1d(res[\"id\"].to_numpy(), target)) == 10\n", +- " \n", +- " ds = ds.create_index(\"vec\", \"IVF_PQ\", metric=metric, num_partitions=num_partitions, num_sub_vectors=num_sub_vectors)\n", +- " \n", ++ "\n", ++ " ds = ds.create_index(\n", ++ " \"vec\",\n", ++ " \"IVF_PQ\",\n", ++ " metric=metric,\n", ++ " num_partitions=num_partitions,\n", ++ " num_sub_vectors=num_sub_vectors,\n", ++ " )\n", ++ "\n", + " recall_data = []\n", + " for nprobes in nprobes_list:\n", + " for refine_factor in refine_factor_list:\n", + " hits = 0\n", + " # check that brute force impl is correct\n", +- " for q, target in zip(tqdm(query, desc=f\"out of sample, nprobes={nprobes}, refine={refine_factor}\"), gt):\n", +- " res = ds.to_table(nearest={\n", +- " \"column\": \"vec\",\n", +- " \"q\": q,\n", +- " \"k\": 10,\n", +- " \"nprobes\": nprobes,\n", +- " \"refine_factor\": refine_factor,\n", +- " }, columns=[\"id\"])[\"id\"].to_numpy()\n", ++ " for q, target in zip(\n", ++ " tqdm(\n", ++ " query,\n", ++ " desc=f\"out of sample, nprobes={nprobes}, refine={refine_factor}\",\n", ++ " ),\n", ++ " gt,\n", ++ " ):\n", ++ " res = ds.to_table(\n", ++ " nearest={\n", ++ " \"column\": \"vec\",\n", ++ " \"q\": q,\n", ++ " \"k\": 10,\n", ++ " \"nprobes\": nprobes,\n", ++ " \"refine_factor\": refine_factor,\n", ++ " },\n", ++ " columns=[\"id\"],\n", ++ " )[\"id\"].to_numpy()\n", + " hits += len(np.intersect1d(res, target))\n", +- " recall_data.append([\n", +- " \"out_of_sample\",\n", +- " nprobes,\n", +- " refine_factor,\n", +- " hits / 10 / len(gt),\n", +- " ])\n", ++ " recall_data.append(\n", ++ " [\n", ++ " \"out_of_sample\",\n", ++ " nprobes,\n", ++ " refine_factor,\n", ++ " hits / 10 / len(gt),\n", ++ " ]\n", ++ " )\n", + " # check that brute force impl is correct\n", +- " for q, target in zip(tqdm(in_sample, desc=f\"in sample nprobes={nprobes}, refine={refine_factor}\"), gt_in_sample):\n", +- " res = ds.to_table(nearest={\n", +- " \"column\": \"vec\",\n", +- " \"q\": q,\n", +- " \"k\": 10,\n", +- " \"nprobes\": nprobes,\n", +- " \"refine_factor\": refine_factor,\n", +- " }, columns=[\"id\"])[\"id\"].to_numpy()\n", ++ " for q, target in zip(\n", ++ " tqdm(\n", ++ " in_sample,\n", ++ " desc=f\"in sample nprobes={nprobes}, refine={refine_factor}\",\n", ++ " ),\n", ++ " gt_in_sample,\n", ++ " ):\n", ++ " res = ds.to_table(\n", ++ " nearest={\n", ++ " \"column\": \"vec\",\n", ++ " \"q\": q,\n", ++ " \"k\": 10,\n", ++ " \"nprobes\": nprobes,\n", ++ " \"refine_factor\": refine_factor,\n", ++ " },\n", ++ " columns=[\"id\"],\n", ++ " )[\"id\"].to_numpy()\n", + " hits += len(np.intersect1d(res, target))\n", +- " recall_data.append([\n", +- " \"in_sample\",\n", +- " nprobes,\n", +- " refine_factor,\n", +- " hits / 10 / len(gt_in_sample),\n", +- " ])\n", ++ " recall_data.append(\n", ++ " [\n", ++ " \"in_sample\",\n", ++ " nprobes,\n", ++ " refine_factor,\n", ++ " hits / 10 / len(gt_in_sample),\n", ++ " ]\n", ++ " )\n", + " return recall_data" + ] + }, +@@ -124,15 +158,19 @@ + "outputs": [], + "source": [ + "def make_plot(recall_data):\n", +- " df = pd.DataFrame(recall_data, columns=[\"case\", \"nprobes\", \"refine_factor\", \"recall\"])\n", +- " \n", ++ " df = pd.DataFrame(\n", ++ " recall_data, columns=[\"case\", \"nprobes\", \"refine_factor\", \"recall\"]\n", ++ " )\n", ++ "\n", + " num_cases = len(df[\"case\"].unique())\n", + " (fig, axs) = plt.subplots(1, 2, figsize=(16, 8))\n", +- " \n", ++ "\n", + " for case, ax in zip(df[\"case\"].unique(), axs):\n", + " current_case = df[df[\"case\"] == case]\n", + " sns.heatmap(\n", +- " current_case.drop(columns=[\"case\"]).set_index([\"nprobes\", \"refine_factor\"])[\"recall\"].unstack(),\n", ++ " current_case.drop(columns=[\"case\"])\n", ++ " .set_index([\"nprobes\", \"refine_factor\"])[\"recall\"]\n", ++ " .unstack(),\n", + " annot=True,\n", + " ax=ax,\n", + " ).set(title=f\"Recall -- {case}\")" diff --git a/test_data/v0.30.0_pre_created_at/datagen.py b/test_data/v0.30.0_pre_created_at/datagen.py index 3d610c1fedc..8a669c12cda 100644 --- a/test_data/v0.30.0_pre_created_at/datagen.py +++ b/test_data/v0.30.0_pre_created_at/datagen.py @@ -2,6 +2,7 @@ Generate test data with Lance 0.29.0 to create an index without created_at field. This tests backward compatibility for the created_at field added in 0.30.0. """ + import lance import pyarrow as pa import pyarrow.compute as pc @@ -13,12 +14,14 @@ ndims = 16 nvecs = 256 -data = pa.table({ - "id": pa.array(range(nvecs)), - "vec": pa.FixedSizeListArray.from_arrays( - pc.random(ndims * nvecs).cast(pa.float32()), ndims - ), -}) +data = pa.table( + { + "id": pa.array(range(nvecs)), + "vec": pa.FixedSizeListArray.from_arrays( + pc.random(ndims * nvecs).cast(pa.float32()), ndims + ), + } +) # Write dataset dataset = lance.write_dataset(data, "index_without_created_at") @@ -33,4 +36,4 @@ print("Created dataset with IVF_PQ index on 'vec' field") print(f"Dataset version: {dataset.version}") -print(f"Index created with Lance {lance.__version__}") \ No newline at end of file +print(f"Index created with Lance {lance.__version__}") diff --git a/test_debug.py b/test_debug.py index e42798d5ff3..3f4452f7bf4 100644 --- a/test_debug.py +++ b/test_debug.py @@ -5,46 +5,52 @@ from lance.namespace import DirectoryNamespace CONFIG = { - 'allow_http': 'true', - 'aws_access_key_id': 'ACCESS_KEY', - 'aws_secret_access_key': 'SECRET_KEY', - 'aws_endpoint': 'http://localhost:4566', - 'aws_region': 'us-east-1', + "allow_http": "true", + "aws_access_key_id": "ACCESS_KEY", + "aws_secret_access_key": "SECRET_KEY", + "aws_endpoint": "http://localhost:4566", + "aws_region": "us-east-1", } storage_options = copy.deepcopy(CONFIG) storage_options_with_refresh = dict(storage_options) -storage_options_with_refresh['refresh_offset_millis'] = '1000' +storage_options_with_refresh["refresh_offset_millis"] = "1000" -dir_props = {f'storage.{k}': v for k, v in storage_options_with_refresh.items()} -dir_props['root'] = 's3://lance-namespace-integtest/namespace_root' -dir_props['ops_metrics_enabled'] = 'true' -dir_props['vend_input_storage_options'] = 'true' -dir_props['vend_input_storage_options_refresh_interval_millis'] = '3600000' +dir_props = {f"storage.{k}": v for k, v in storage_options_with_refresh.items()} +dir_props["root"] = "s3://lance-namespace-integtest/namespace_root" +dir_props["ops_metrics_enabled"] = "true" +dir_props["vend_input_storage_options"] = "true" +dir_props["vend_input_storage_options_refresh_interval_millis"] = "3600000" namespace = DirectoryNamespace(**dir_props) -table1 = pa.Table.from_pylist([{'a': 1, 'b': 2}]) -table_name = 'debug_print3_' + uuid.uuid4().hex -table_id = ['test_ns', table_name] +table1 = pa.Table.from_pylist([{"a": 1, "b": 2}]) +table_name = "debug_print3_" + uuid.uuid4().hex +table_id = ["test_ns", table_name] -print('=== Creating table ===') +print("=== Creating table ===") ds = lance.write_dataset( table1, namespace_client=namespace, table_id=table_id, - mode='create', + mode="create", storage_options=storage_options, ) -print('=== Table created ===') -print('Describe count after write:', namespace.retrieve_ops_metrics().get('describe_table', 0)) +print("=== Table created ===") +print( + "Describe count after write:", + namespace.retrieve_ops_metrics().get("describe_table", 0), +) -print('') -print('=== Opening dataset via lance.dataset() ===') +print("") +print("=== Opening dataset via lance.dataset() ===") ds_from_namespace = lance.dataset( namespace_client=namespace, table_id=table_id, storage_options=storage_options, ) -print('=== Dataset opened ===') -print('Describe count after lance.dataset():', namespace.retrieve_ops_metrics().get('describe_table', 0)) +print("=== Dataset opened ===") +print( + "Describe count after lance.dataset():", + namespace.retrieve_ops_metrics().get("describe_table", 0), +) From b54450b8fb3f9893d3d7482d7dfc71c9f537a766 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 18:57:11 -0400 Subject: [PATCH 03/31] =?UTF-8?q?Set=20=CF=83=E1=B5=A2=20=3D=20=E2=88=9A?= =?UTF-8?q?=CE=BB=E1=B5=A2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- rust/lance-linalg/src/svd.rs | 171 +---------------------------------- 1 file changed, 1 insertion(+), 170 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 66732ffcd3e..676a14d1dc9 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -1,25 +1,6 @@ const EPS: f64 = 2.220_446_049_250_313e-16; const MAX_JACOBI_SWEEPS: usize = 200; -/// x² — just multiplication, listed for clarity -#[inline] -fn sq(x: f64) -> f64 { x * x } - -/// Raise f64 to an integer power. -fn powi(mut x: f64, mut n: i32) -> f64 { - if n == 0 { return 1.0; } - if n < 0 { x = 1.0 / x; n = -n; } - let mut result = 1.0f64; - let mut base = x; - let mut exp = n as u32; - while exp > 0 { - if exp & 1 == 1 { result *= base; } - base *= base; - exp >>= 1; - } - result -} - fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { let mut s = a.to_vec(); let mut v = eye(n); @@ -59,80 +40,6 @@ fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { (eigenvalues, v) } -fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { - let s_pp = s[p * n + p]; - let s_qq = s[q * n + q]; - let s_pq = s[p * n + q]; - - s[p * n + p] = sq(c) * s_pp - 2.0 * sn * c * s_pq + sq(sn) * s_qq; - s[q * n + q] = sq(sn) * s_pp + 2.0 * sn * c * s_pq + sq(c) * s_qq; - s[p * n + q] = 0.0; - s[q * n + p] = 0.0; - - for r in 0..n { - if r == p || r == q { continue; } - let s_rp = s[r * n + p]; - let s_rq = s[r * n + q]; - let new_rp = c * s_rp - sn * s_rq; - let new_rq = sn * s_rp + c * s_rq; - s[r * n + p] = new_rp; s[p * n + r] = new_rp; - s[r * n + q] = new_rq; s[q * n + r] = new_rq; - } -} - -fn apply_givens_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { - for r in 0..n { - let vp = v[r * n + p]; - let vq = v[r * n + q]; - v[r * n + p] = c * vp - sn * vq; - v[r * n + q] = sn * vp + c * vq; - } -} - -// ============================================================ -// Modified Gram-Schmidt orthonormalization -// ============================================================ - -fn gram_schmidt(cols: &mut Vec>) { - let ncols = cols.len(); - for i in 0..ncols { - for j in 0..i { - let dot: f64 = cols[i].iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); - let cj = cols[j].clone(); - for (a, b) in cols[i].iter_mut().zip(cj.iter()) { - *a -= dot * b; - } - } - let norm = (cols[i].iter().map(|&x| sq(x)).sum::()).sqrt(); - if norm > EPS { - for x in cols[i].iter_mut() { *x /= norm; } - } else { - // Replace with an orthonormal basis vector not already spanned - let dim = cols[i].len(); - 'search: for k in 0..dim { - let mut e = vec![0.0f64; dim]; - e[k] = 1.0; - for j in 0..i { - let dot: f64 = e.iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); - let cj = cols[j].clone(); - for (a, b) in e.iter_mut().zip(cj.iter()) { *a -= dot * b; } - } - let n2 = (e.iter().map(|&x| sq(x)).sum::()).sqrt(); - if n2 > EPS { - for x in e.iter_mut() { *x /= n2; } - cols[i] = e; - break 'search; - } - } - } - } -} - -// ============================================================ -// Matrix helpers -// ============================================================ - -/// C = Aᵀ·A where A is m×n → C is n×n fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { let mut c = vec![0f64; n * n]; for i in 0..n { @@ -146,90 +53,14 @@ fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { c } -/// y = A·x (A is m×n row-major, x length n → y length m) -fn mat_vec_mul(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { - (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() -} - -fn eye(n: usize) -> Vec { - let mut m = vec![0f64; n * n]; - for i in 0..n { m[i * n + i] = 1.0; } - m -} - pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { let ata = mat_mul_atb(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); - + let k = m.min(n); let mut sigma: Vec = eigenvalues .iter() .map(|&l| (if l > 0.0 { l } else { 0.0 }).sqrt()) .collect(); - let mut u_cols: Vec> = Vec::with_capacity(m); - for i in 0..k { - let vi: Vec = (0..n).map(|r| v[r * n + i]).collect(); - let av = mat_vec_mul(a, &vi, m, n); - if sigma[i] > EPS * 10.0 { - u_cols.push(av.iter().map(|&x| x / sigma[i]).collect()); - } else { - sigma[i] = 0.0; - u_cols.push(vec![0.0; m]); - } - } - for _ in k..m { - u_cols.push(vec![0.0; m]); - } - gram_schmidt(&mut u_cols); - - let mut u = vec![0f64; m * m]; - for (ci, col) in u_cols.iter().enumerate() { - for (ri, &value) in col.iter().enumerate() { - u[ri * m + ci] = value; // was: ri * m — m was undefined - } - } - - let mut order: Vec = (0..k).collect(); - for i in 1..k { - let mut j = i; - while j > 0 && sigma[order[j - 1]] < sigma[order[j]] { - order.swap(j - 1, j); - j -= 1; - } - } - - let sigma_sorted: Vec = order.iter().map(|&i| sigma[i]).collect(); // was: .map[( - - let mut u_sorted = vec![0f64; m * m]; - for (new_col, &old_col) in order.iter().enumerate() { - for r in 0..m { - u_sorted[r * m + new_col] = u[r * m + old_col]; // was: r * m (m undefined) - } - } - for col in k..m { - for r in 0..m { - u_sorted[r * m + col] = u[r * m + col]; // was: k..m (m undefined) - } - } - - let mut vt = vec![0f64; n * n]; - for (new_row, &old_col) in order.iter().enumerate() { - for c in 0..n { - vt[new_row * n + c] = v[c * n + old_col]; // was: V_transpose, n undefined - } - } - for row in k..n { - for c in 0..n { - vt[row * n + c] = v[c * n + row]; // was: n undefined - } - } - - (u_sorted, sigma_sorted, vt) -} - -fn main() { - let a = vec![3.0, 2.0, 2.0, 2.0, 3.0, -2.0]; - let (u, s, vt) = svd(&a, 2, 3); - println!("{:?}", s); } \ No newline at end of file From 0cea39f1b3ca714f20f1042ac0d88c67688abc02 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 19:20:42 -0400 Subject: [PATCH 04/31] =?UTF-8?q?Sorted=20all=20n=20eigenvalues=20in=20des?= =?UTF-8?q?cending=20order=20before=20computing=20singular=20values=20?= =?UTF-8?q?=CF=83=E1=B5=A2=20=3D=20=E2=88=9A=CE=BB=E1=B5=A2=20for=20the=20?= =?UTF-8?q?top=20k=20eigenvalues?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- rust/lance-linalg/src/svd.rs | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 676a14d1dc9..088e118a627 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -58,9 +58,19 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { let (eigenvalues, v) = jacobi_eigen(&ata, n); let k = m.min(n); - let mut sigma: Vec = eigenvalues + + let mut order: Vec = (0..n).collect(); + for i in 1..n { + let mut j = i; + while j > 0 && eigenvalues[order[j - 1]] < eigenvalues[order[j]] { + order.swap(j - 1, j); + j -= 1; + } + } + + let mut sigma: Vec = order[..k] .iter() - .map(|&l| (if l > 0.0 { l } else { 0.0 }).sqrt()) + .map(|&i| if eigenvalues[i] > 0.0 { eigenvalues[i].sqrt() } else { 0.0 }) .collect(); } \ No newline at end of file From e26009ac6cd75d77264131ca18b71f94e3decb62 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 19:36:18 -0400 Subject: [PATCH 05/31] =?UTF-8?q?Computed=20left=20singular=20vectors=20u?= =?UTF-8?q?=E1=B5=A2=20=3D=20(1/=CF=83=E1=B5=A2)=20*=20A=20*=20v=E1=B5=A2.?= =?UTF-8?q?=20Used=20the=20sorted=20eigenvectors=20to=20build=20the=20U=20?= =?UTF-8?q?matrix.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- rust/lance-linalg/src/svd.rs | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 088e118a627..b6ae56b389d 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -53,6 +53,10 @@ fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { c } +fn mat_vec_mul(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { + (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect(); +} + pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { let ata = mat_mul_atb(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); @@ -73,4 +77,16 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { .map(|&i| if eigenvalues[i] > 0.0 { eigenvalues[i].sqrt() } else { 0.0 }) .collect(); + let mut u_cols: Vec> = Vec::with_capacity(m); + for index in 0..k { + let ei = order[index]; + let vi: Vec = (0..n).map(|r| v[r * n + ei]).collect(); + let av = mat_vec_mul(a, &vi, m, n); + if sigma[index] > EPS * 10.0 { + u_cols.push(av.iter().map(|&x| x / sigma[index]).collect()); + } else { + sigma[index] = 0.0; + u_cols.push(vec![0.0; m]); + } + } } \ No newline at end of file From 71fc8cb55cac34382cfada82e90aaf6088bafede Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 20:11:56 -0400 Subject: [PATCH 06/31] Padded the U matrix with zero columns if m > k. Used Gram-Schmidt to fill these zero columns with orthonormal vectors that complete the basis of R^m. --- rust/lance-linalg/src/svd.rs | 46 ++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index b6ae56b389d..3985d3258e4 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -53,6 +53,47 @@ fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { c } +fn gram_schmidt(cols: &mut Vec>) { + let n_cols = cols.len(); + for i in 0..ncols { + for j in 0..i { + let dot: f64 = cols[i].iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); + let cj = cols[j].clone(); + for (a, b) in cols[i].iter_mut().zip(cj.iter()) { + *a -= dot * b; + } + } + + let norm = cols[i].iter().map(|&x| x * x).sum::().sqrt(); + if norm > EPS { + for x in cols[i].iter_mut() { + *x /= norm; + } + } else { + let dim = cols[i].len(); + 'search: for k in 0..dim { + let mut e = vec![0.0f64; dim]; + e[k] = 1.0; + for j in 0..i { + let dot: f64 = e.iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); + let cj = cols[j].clone(); + for (a, b) in e.iter_mut().zip(cj.iter()) { + *a -= dot * b; + } + } + let n2 = e.iter().map(|&x| x * x).sum::().sqrt(); + if n2 > EPS { + for x in e.iter_mut() { + *x /= n2; + } + cols[i] = e; + break 'search; + } + } + } + } +} + fn mat_vec_mul(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect(); } @@ -89,4 +130,9 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { u_cols.push(vec![0.0; m]); } } + + for _ in k..m { + u_cols.push(vec![0.0; m]); + } + gram_schmidt(&mut u_cols); } \ No newline at end of file From 812d848e0c5b4a7114023bfa961de6faa9383d5d Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 21:09:38 -0400 Subject: [PATCH 07/31] Packed the columns of the U matrix into a row-major mxm matrix --- rust/lance-linalg/src/svd.rs | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 3985d3258e4..fb5880f60b9 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -135,4 +135,11 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { u_cols.push(vec![0.0; m]); } gram_schmidt(&mut u_cols); + + let mut u = vec![0f64; m * m]; + for (ci, col) in u_cols.iter().enumerate() { + for (ri, &value) in col.iter().enumerate() { + u[ri * m + ci] = value; + } + } } \ No newline at end of file From 52972c8d68e279544daac006bbc435de4e746a6e Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 21:25:28 -0400 Subject: [PATCH 08/31] Built the V^T matrix and returned the U, sigma, and V^T matrices --- rust/lance-linalg/src/svd.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index fb5880f60b9..75eaa87941b 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -142,4 +142,14 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { u[ri * m + ci] = value; } } + + let mut vt = vec![0f64; n * n]; + for new_row in 0..n { + let old_col = order[new_row]; + for c in 0..n { + vt[new_row * n + c] = v[c * n + old_col]; + } + } + + (u, sigma, vt) } \ No newline at end of file From 2b347fe622fc39e70124cdb96b1b5aed4eaa0750 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 21:40:54 -0400 Subject: [PATCH 09/31] Implemented the create_identity_matrix() function, which creates an nxn identity matrix as a row-major vector --- rust/lance-linalg/src/svd.rs | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 75eaa87941b..7e3a0fce42c 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -1,9 +1,25 @@ const EPS: f64 = 2.220_446_049_250_313e-16; const MAX_JACOBI_SWEEPS: usize = 200; +fn create_identity_matrix(n: usize) -> Vec { + let mut m = vec![0f64; n * n]; + for i in 0..n { + m[i * n + i] = 1.0; + } + m +} + +fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + +} + +fn apply_givens_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + +} + fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { let mut s = a.to_vec(); - let mut v = eye(n); + let mut v = create_identity_matrix(n); for _ in 0..MAX_JACOBI_SWEEPS { let mut max_val = 0.0f64; From 7202b70a79364a5240a01e95bb46a3a7a1f9bdfd Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 22:58:41 -0400 Subject: [PATCH 10/31] Created the jacobi_rotate() function, which applies a Givens rotation to a symmetric matrix S. Also created the apply_givens_right() function, which accumulates a Givens rotation into the matrix V from the right --- rust/lance-linalg/src/svd.rs | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 7e3a0fce42c..8663852b3bc 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -10,11 +10,39 @@ fn create_identity_matrix(n: usize) -> Vec { } fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { - + let s_pp = s[p * n + p]; + let s_qq = s[q * n + q]; + let s_pq = s[p * n + q]; + + s[p * n + p] = c * c * s_pp - 2.0 * sn * c * s_pq + sn * sn * s_qq; + s[q * n + q] = sn * sn * s_pp + 2.0 * sn * c * s_pq + c * c * s_qq; + s[p * n + q] = 0.0; + s[q * n + p] = 0.0; + + for r in 0..n { + if r == p || r == q { + continue; + } + let s_rp = s[r * n + p]; + let s_rq = s[r * n + q]; + let new_rp = c * s_rp - sn * s_rq; + let new_rq = sn * s_rp + c * s_rq; + + s[r * n + p] = new_rp; + s[p * n + r] = new_rp; + s[r * n + q] = new_rq; + s[q * n + r] = new_rq; + } } fn apply_givens_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + for r in 0..n { + let vp = v[r * n + p]; + let vq = v[r * n + q]; + v[r * n + p] = c * vp - sn * vq; + v[r * n + q] = sn * vp + c * vq; + } } fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { From 49c6009ccac477350269a190210989006b436319 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:02:18 -0400 Subject: [PATCH 11/31] Renamed a few functions for clarity purposes --- rust/lance-linalg/src/svd.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 8663852b3bc..50ba3705dd2 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -35,7 +35,7 @@ fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { } } -fn apply_givens_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { +fn apply_givens_rotation_from_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { for r in 0..n { let vp = v[r * n + p]; let vq = v[r * n + q]; @@ -77,14 +77,14 @@ fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { let (sn, c) = theta.sin_cos(); jacobi_rotate(&mut s, n, p, q, c, sn); - apply_givens_right(&mut v, n, p, q, c, sn); + apply_givens_rotation_from_right(&mut v, n, p, q, c, sn); } let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); (eigenvalues, v) } -fn mat_mul_atb(a: &[f64], m: usize, n: usize) -> Vec { +fn compute_ata(a: &[f64], m: usize, n: usize) -> Vec { let mut c = vec![0f64; n * n]; for i in 0..n { for l in 0..m { @@ -138,12 +138,12 @@ fn gram_schmidt(cols: &mut Vec>) { } } -fn mat_vec_mul(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { +fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect(); } pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { - let ata = mat_mul_atb(a, m, n); + let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); let k = m.min(n); @@ -166,7 +166,7 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { for index in 0..k { let ei = order[index]; let vi: Vec = (0..n).map(|r| v[r * n + ei]).collect(); - let av = mat_vec_mul(a, &vi, m, n); + let av = multiply_A_by_vector(a, &vi, m, n); if sigma[index] > EPS * 10.0 { u_cols.push(av.iter().map(|&x| x / sigma[index]).collect()); } else { From 3d64e1f036baa5b711f51c5a392eb0c65cb4d896 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:05:51 -0400 Subject: [PATCH 12/31] Fixed a typo in a variable name --- rust/lance-linalg/src/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 50ba3705dd2..960cd3801d5 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -98,7 +98,7 @@ fn compute_ata(a: &[f64], m: usize, n: usize) -> Vec { } fn gram_schmidt(cols: &mut Vec>) { - let n_cols = cols.len(); + let ncols = cols.len(); for i in 0..ncols { for j in 0..i { let dot: f64 = cols[i].iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); From 5602eedd5179f59912ae0b691acd82dd2368ea6c Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:09:07 -0400 Subject: [PATCH 13/31] Removed a semicolon from a return value --- rust/lance-linalg/src/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 960cd3801d5..11ee13bb8f6 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -139,7 +139,7 @@ fn gram_schmidt(cols: &mut Vec>) { } fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { - (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect(); + (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() } pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { From 4dd41f8959b14b6ec9fa7da57cf17577661cc6e9 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:44:16 -0400 Subject: [PATCH 14/31] Checked whether the input matrix has at least one row and at least one column --- rust/lance-linalg/src/svd.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 11ee13bb8f6..ad11e7c84ed 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -143,6 +143,8 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { } pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { + assert!(m > 0 && n > 0, "Matrix must have at least 1 row and at least 1 column."); + let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); From 74fed31a20bb1ce88330e3b5e2714b2dc7075031 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:45:43 -0400 Subject: [PATCH 15/31] Checked whether the data length of the input matrix matches the product of the specified number of rows and number of columns --- rust/lance-linalg/src/svd.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index ad11e7c84ed..6f68580dfc4 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -144,6 +144,7 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { assert!(m > 0 && n > 0, "Matrix must have at least 1 row and at least 1 column."); + assert_eq!(a.len(), m * n, "Data length of matrix does not match the product of the specified number of rows and number of columns", a.len(), m, n); let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); From 80fd348ff5f652594a2bb2fd60e12612efbc277b Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 22 Jun 2026 23:47:28 -0400 Subject: [PATCH 16/31] Checked whether the input matrix contains null or infinite entries --- rust/lance-linalg/src/svd.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 6f68580dfc4..5dabcf7eb57 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -145,6 +145,7 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { assert!(m > 0 && n > 0, "Matrix must have at least 1 row and at least 1 column."); assert_eq!(a.len(), m * n, "Data length of matrix does not match the product of the specified number of rows and number of columns", a.len(), m, n); + assert!(a.iter().all(|x| x.is_finite()), "Matrix contains null or infinite entries."); let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); From 6187fde1f6a78212e0305573a3fd9e0c0273af9c Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 00:01:25 -0400 Subject: [PATCH 17/31] Edited my input validation statements for clarity purposes --- rust/lance-linalg/src/svd.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 5dabcf7eb57..642fd86d81b 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -144,8 +144,8 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { assert!(m > 0 && n > 0, "Matrix must have at least 1 row and at least 1 column."); - assert_eq!(a.len(), m * n, "Data length of matrix does not match the product of the specified number of rows and number of columns", a.len(), m, n); - assert!(a.iter().all(|x| x.is_finite()), "Matrix contains null or infinite entries."); + assert_eq!(a.len(), m * n, "Data length of matrix must match the product of the specified number of rows and number of columns.", a.len(), m, n); + assert!(a.iter().all(|x| x.is_finite()), "Matrix must not contain null or infinite entries."); let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); From ff2e2720b79df096c00411e613c12248e64405b0 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 17:44:33 -0400 Subject: [PATCH 18/31] Converted my input validation statement for whether the matrix has at least one row and at least one column, from an assert statement to an if statement. --- rust/lance-linalg/src/svd.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 642fd86d81b..c7588d7077f 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -143,7 +143,10 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { } pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { - assert!(m > 0 && n > 0, "Matrix must have at least 1 row and at least 1 column."); + if m == 0 || n == 0 { + println!("Error: Matrix must have at least 1 row and at least 1 column."); + return (vec![], vec![], vec![]); + } assert_eq!(a.len(), m * n, "Data length of matrix must match the product of the specified number of rows and number of columns.", a.len(), m, n); assert!(a.iter().all(|x| x.is_finite()), "Matrix must not contain null or infinite entries."); From 7bb18c1d6c71a4d478383ef8e3a4ad9c38619376 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 17:46:51 -0400 Subject: [PATCH 19/31] Converted my input validation statement for whether the data length of the matrix matches the product of the specified number of rows and number of columns, from an assert statement to an if statement. --- rust/lance-linalg/src/svd.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index c7588d7077f..001d54aafd4 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -147,7 +147,10 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { println!("Error: Matrix must have at least 1 row and at least 1 column."); return (vec![], vec![], vec![]); } - assert_eq!(a.len(), m * n, "Data length of matrix must match the product of the specified number of rows and number of columns.", a.len(), m, n); + if a.len() != m * n { + println!("Data length of matrix must match the product of the specified number of rows and number of columns."); + return (vec![], vec![], vec![]); + } assert!(a.iter().all(|x| x.is_finite()), "Matrix must not contain null or infinite entries."); let ata = compute_ata(a, m, n); From 93584c82be75a9cfd2098faf09ead70ee353f2b7 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 17:48:45 -0400 Subject: [PATCH 20/31] Converted my input validation statement for whether the matrix contains null or infinite entries, from an assert statement to an if statement. --- rust/lance-linalg/src/svd.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 001d54aafd4..4e8a5935f0d 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -151,7 +151,10 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { println!("Data length of matrix must match the product of the specified number of rows and number of columns."); return (vec![], vec![], vec![]); } - assert!(a.iter().all(|x| x.is_finite()), "Matrix must not contain null or infinite entries."); + if !a.iter().all(|x| x.is_finite()) { + println!("Matrix must not contain null or infinite entries."); + return (vec![], vec![], vec![]); + } let ata = compute_ata(a, m, n); let (eigenvalues, v) = jacobi_eigen(&ata, n); From fe78e53490fb4f7413755ccc9d132f1665960803 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 17:50:31 -0400 Subject: [PATCH 21/31] Edited all my input validation logic so that multiple different error messages can be displayed if needed --- rust/lance-linalg/src/svd.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 4e8a5935f0d..b15100ae868 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -145,14 +145,15 @@ fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { if m == 0 || n == 0 { println!("Error: Matrix must have at least 1 row and at least 1 column."); - return (vec![], vec![], vec![]); } if a.len() != m * n { println!("Data length of matrix must match the product of the specified number of rows and number of columns."); - return (vec![], vec![], vec![]); } if !a.iter().all(|x| x.is_finite()) { println!("Matrix must not contain null or infinite entries."); + } + + if m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { return (vec![], vec![], vec![]); } From bcfca5cf47638e170363bc0de08c261a0098f785 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Tue, 23 Jun 2026 20:39:24 -0400 Subject: [PATCH 22/31] Included unit, integration, and manual tests and made a minor edit to my input validation error messages --- rust/lance-linalg/src/svd.rs | 292 ++++++++++++++++++++++++++++++++++- 1 file changed, 288 insertions(+), 4 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index b15100ae868..800a00f0ac8 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -138,7 +138,7 @@ fn gram_schmidt(cols: &mut Vec>) { } } -fn multiply_A_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { +fn multiply_a_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() } @@ -147,10 +147,10 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { println!("Error: Matrix must have at least 1 row and at least 1 column."); } if a.len() != m * n { - println!("Data length of matrix must match the product of the specified number of rows and number of columns."); + println!("Error: Data length of matrix must match the product of the specified number of rows and number of columns."); } if !a.iter().all(|x| x.is_finite()) { - println!("Matrix must not contain null or infinite entries."); + println!("Error: Matrix must not contain null or infinite entries."); } if m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { @@ -180,7 +180,7 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { for index in 0..k { let ei = order[index]; let vi: Vec = (0..n).map(|r| v[r * n + ei]).collect(); - let av = multiply_A_by_vector(a, &vi, m, n); + let av = multiply_a_by_vector(a, &vi, m, n); if sigma[index] > EPS * 10.0 { u_cols.push(av.iter().map(|&x| x / sigma[index]).collect()); } else { @@ -210,4 +210,288 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { } (u, sigma, vt) +} + +fn main() { + println!("========== UNIT TESTS =========="); + unit_test_1_multiply_a_by_vector(); + unit_test_2_compute_ata(); + unit_test_3_create_identity_matrix(); + + println!("\n========== INTEGRATION TESTS =========="); + integration_test_1_svd_2x3(); + integration_test_2_svd_4x4(); + + println!("\n========== MANUAL TESTS =========="); + manual_test_1_svd_3x3(); + println!(); + manual_test_2_svd_empty_matrix(); + println!(); + manual_test_3_svd_length_mismatch(); + println!(); + manual_test_4_svd_nan_entry(); + println!(); + manual_test_5_svd_infinite_entry(); + println!(); + manual_test_6_jacobi_rotate(); + println!(); + manual_test_7_apply_givens_rotation_from_right(); + println!(); + manual_test_8_jacobi_eigen(); + println!(); + manual_test_9_gram_schmidt(); + + println!("\nAll tests complete."); +} + +// ========================================================================= +// UNIT TESTS — automated +// ========================================================================= + +fn unit_test_1_multiply_a_by_vector() { + let a = vec![5.0, 7.0, -1.0, 4.0]; + let x = vec![-9.0, -2.0]; + let m = 2; + let n = 2; + let result = multiply_a_by_vector(&a, &x, m, n); + assert_eq!(result, vec![-59.0, 1.0], + "unit_test_1 FAILED: got {:?}, expected [-59.0, 1.0]", result); + println!("unit_test_1 PASSED"); +} + +fn unit_test_2_compute_ata() { + let a = vec![3.0, 2.0, -4.0, 7.0, -2.0, -1.0, 5.0, 4.0]; + let m = 2; + let n = 4; + let result = compute_ata(&a, m, n); + let expected = vec![ + 13.0, 8.0, -22.0, 13.0, + 8.0, 5.0, -13.0, 10.0, + -22.0, -13.0, 41.0, -8.0, + 13.0, 10.0, -8.0, 65.0, + ]; + for (i, (r, e)) in result.iter().zip(expected.iter()).enumerate() { + assert!( + (r - e).abs() < 1e-9, + "unit_test_2 FAILED at index {}: got {}, expected {}", i, r, e + ); + } + println!("unit_test_2 PASSED"); +} + +fn unit_test_3_create_identity_matrix() { + let n = 5; + let result = create_identity_matrix(n); + let expected = vec![ + 1.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0, + ]; + assert_eq!(result, expected, + "unit_test_3 FAILED: got {:?}", result); + println!("unit_test_3 PASSED"); +} + +// ========================================================================= +// INTEGRATION TESTS — automated +// ========================================================================= + +fn integration_test_1_svd_2x3() { + let a = vec![3.0, 2.0, 2.0, 2.0, 3.0, -2.0]; + let m = 2; + let n = 3; + let (u, sigma, vt) = svd(&a, m, n); + + assert!(!sigma.is_empty(), "integration_test_1 FAILED: sigma is empty"); + assert!( + (sigma[0] - 5.0).abs() < 1e-3, + "integration_test_1 FAILED: sigma[0] = {}, expected 5.0", sigma[0] + ); + assert!( + (sigma[1] - 3.0).abs() < 1e-3, + "integration_test_1 FAILED: sigma[1] = {}, expected 3.0", sigma[1] + ); + + // : Vec added — this is the only change, asserts are unchanged + let expected_u: Vec = vec![0.7071, 0.7071, 0.7071, -0.7071]; + let expected_vt: Vec = vec![ + 0.7071, 0.7071, 0.0000, + 0.2357, -0.2357, 0.9428, + 0.6667, -0.6667, -0.3333, + ]; + for (i, (&r, &e)) in u.iter().zip(expected_u.iter()).enumerate() { + assert!( + (r.abs() - e.abs()).abs() < 1e-3, + "integration_test_1 FAILED U at index {}: got {}, expected {}", i, r, e + ); + } + for (i, (&r, &e)) in vt.iter().zip(expected_vt.iter()).enumerate() { + assert!( + (r.abs() - e.abs()).abs() < 1e-3, + "integration_test_1 FAILED V^T at index {}: got {}, expected {}", i, r, e + ); + } + println!("integration_test_1 PASSED: U = {:?}, sigma = {:?}, V^T = {:?}", u, sigma, vt); +} + +fn integration_test_2_svd_4x4() { + let a = vec![ + 7.0, -4.0, 5.0, 5.0, + 8.0, -2.0,-10.0, 1.0, + -1.0, -8.0, 9.0, 3.0, + 8.0, 7.0, -3.0, 4.0, + ]; + let m = 4; + let n = 4; + let (u, sigma, vt) = svd(&a, m, n); + + assert!(!sigma.is_empty(), "integration_test_2 FAILED: sigma is empty"); + + // : Vec added — asserts unchanged + let expected_sigma: Vec = vec![17.834861, 13.682368, 8.433152, 0.769723]; + for (i, (&r, &e)) in sigma.iter().zip(expected_sigma.iter()).enumerate() { + assert!( + (r - e).abs() < 1e-3, + "integration_test_2 FAILED sigma at index {}: got {}, expected {}", i, r, e + ); + } + + let expected_u: Vec = vec![ + 0.1103, 0.7669, 0.1056, -0.6233, + -0.5973, 0.3163, -0.7191, 0.1617, + 0.5994, 0.4553, -0.1537, 0.6402, + -0.5214, 0.3233, 0.6695, 0.4189, + ]; + let expected_vt: Vec = vec![ + -0.4921, -0.4312, 0.7560, -0.0187, + 0.7330, -0.3713, 0.2777, 0.4977, + 0.0588, 0.8219, 0.5131, 0.2402, + -0.4659, -0.0249, -0.2968, 0.8332, + ]; + for (i, (&r, &e)) in u.iter().zip(expected_u.iter()).enumerate() { + assert!( + (r.abs() - e.abs()).abs() < 1e-3, + "integration_test_2 FAILED U at index {}: got {}, expected {}", i, r, e + ); + } + for (i, (&r, &e)) in vt.iter().zip(expected_vt.iter()).enumerate() { + assert!( + (r.abs() - e.abs()).abs() < 1e-3, + "integration_test_2 FAILED V^T at index {}: got {}, expected {}", i, r, e + ); + } + println!("integration_test_1 PASSED: U = {:?}, sigma = {:?}, V^T = {:?}", u, sigma, vt); +} + +// ========================================================================= +// MANUAL TESTS +// ========================================================================= + +fn manual_test_1_svd_3x3() { + let a = vec![-9.0, -5.0, -2.0, 4.0, -1.0, 6.0, 9.0, -2.0, -6.0]; + let m = 3; + let n = 3; + let (u, sigma, vt) = svd(&a, m, n); + println!("manual_test_1 U: {:?}", u); + println!("manual_test_1 sigma: {:?}", sigma); + println!("manual_test_1 V^T: {:?}", vt); + println!("Expected U ≈ [-0.6955, -0.4291, -0.5763, 0.2417, 0.6157, -0.7500, 0.6767, -0.6609, -0.3245], Expected sigma ≈ [13.5037, 8.9903, 4.5633], Expected V^T ≈ [0.9861, 0.1394, -0.0903, 0.0419, 0.3172, 0.9474, -0.1607, 0.9381, -0.3070]"); +} + +fn manual_test_2_svd_empty_matrix() { + let a: Vec = vec![]; + let m = 0; + let n = 0; + let (u, sigma, vt) = svd(&a, m, n); + println!("manual_test_2 U: {:?}", u); + println!("manual_test_2 sigma: {:?}", sigma); + println!("manual_test_2 V^T: {:?}", vt); + println!("Expected: all empty vectors, error message: 'Error: Matrix must have at least 1 row and at least 1 column.'"); +} + +fn manual_test_3_svd_length_mismatch() { + let a = vec![9.0, -8.0, 3.0, -1.0]; + let m = 1; + let n = 2; + let (u, sigma, vt) = svd(&a, m, n); + println!("manual_test_3 sigma: {:?}", sigma); + println!("manual_test_3 U: {:?}", u); + println!("manual_test_3 V^T: {:?}", vt); + println!("Expected: all empty vectors, error message: 'Error: Data length of matrix must match the product of the specified number of rows and number of columns.'"); +} + +fn manual_test_4_svd_nan_entry() { + let a = vec![3.0, 1.0, f64::NAN, -4.0, -2.0, -1.0, 8.0, 3.0, 1.0]; + let m = 3; + let n = 3; + let (u, sigma, vt) = svd(&a, m, n); + println!("manual_test_4 sigma: {:?}", sigma); + println!("manual_test_4 U: {:?}", u); + println!("manual_test_4 V^T: {:?}", vt); + println!("Expected: all empty vectors, error message: 'Error: Matrix must not contain null or infinite entries.'"); +} + +fn manual_test_5_svd_infinite_entry() { + let a = vec![9.0, f64::INFINITY, -1.0, 8.0]; + let m = 2; + let n = 2; + let (u, sigma, vt) = svd(&a, m, n); + println!("manual_test_5 sigma: {:?}", sigma); + println!("manual_test_5 U: {:?}", u); + println!("manual_test_5 V^T: {:?}", vt); + println!("Expected: all empty vectors, error message: 'Error: Matrix must not contain null or infinite entries.'"); +} + +fn manual_test_6_jacobi_rotate() { + let mut a = vec![4.0, 2.0, 2.0, 3.0]; + let p = 0; + let q = 1; + let s_pq: f64 = a[p * 2 + q]; + let diff: f64 = a[q * 2 + q] - a[p * 2 + p]; + let theta: f64 = 0.5 * (2.0 * s_pq / diff).atan(); + let (sn, c) = theta.sin_cos(); + jacobi_rotate(&mut a, 2, p, q, c, sn); + println!("manual_test_6 after jacobi_rotate: {:?}", a); + println!("Off-diagonal a[0,1] = {:.6} (expected ~0.0)", a[p * 2 + q]); + println!("Trace = {:.6} (expected 7.0)", a[0] + a[3]); +} + +fn manual_test_7_apply_givens_rotation_from_right() { + let mut b = vec![1.0, 0.0, 0.0, 1.0]; + let angle: f64 = std::f64::consts::FRAC_PI_4; + let c = angle.cos(); + let s = angle.sin(); + apply_givens_rotation_from_right(&mut b, 2, 0, 1, c, s); + println!("manual_test_7 after rotation: {:?}", b); + println!("Expected: [0.7071, -0.7071, 0.7071, 0.7071]"); +} + +fn manual_test_8_jacobi_eigen() { + let a = vec![4.0, 1.0, 1.0, 3.0]; + let (eigenvalues, eigenvectors) = jacobi_eigen(&a, 2); + println!("manual_test_8 eigenvalues: {:?}", eigenvalues); + println!("manual_test_8 eigenvectors: {:?}", eigenvectors); + println!("Expected eigenvalues ≈ [4.618, 2.382] (any order)"); +} + +fn manual_test_9_gram_schmidt() { + let mut cols = vec![ + vec![1.0, 1.0, 0.0], + vec![1.0, 0.0, 1.0], + vec![0.0, 1.0, 1.0], + ]; + gram_schmidt(&mut cols); + println!("manual_test_9 after gram_schmidt:"); + for (i, col) in cols.iter().enumerate() { + let norm: f64 = col.iter().map(|&x| x * x).sum::().sqrt(); + println!(" col[{}] = {:?} norm = {:.6}", i, col, norm); + } + let dot01: f64 = cols[0].iter().zip(cols[1].iter()).map(|(&a, &b)| a * b).sum(); + let dot02: f64 = cols[0].iter().zip(cols[2].iter()).map(|(&a, &b)| a * b).sum(); + let dot12: f64 = cols[1].iter().zip(cols[2].iter()).map(|(&a, &b)| a * b).sum(); + println!("dot(col0,col1) = {:.6} (expected ~0.0)", dot01); + println!("dot(col0,col2) = {:.6} (expected ~0.0)", dot02); + println!("dot(col1,col2) = {:.6} (expected ~0.0)", dot12); } \ No newline at end of file From 355603b12cc7168f5ec4d33afb00753a02ee7b3c Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 29 Jun 2026 23:01:19 -0400 Subject: [PATCH 23/31] Edited my input validation logic so that the error message on empty input matrices will be displayed if the length of the input matrix is 0, even if the number of specified rows and columns is not 0. --- rust/lance-linalg/src/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 800a00f0ac8..c2b09dc6ebf 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -143,7 +143,7 @@ fn multiply_a_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { } pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { - if m == 0 || n == 0 { + if a.len() == 0 || m == 0 || n == 0 { println!("Error: Matrix must have at least 1 row and at least 1 column."); } if a.len() != m * n { From 70b006ccd8b7f250fbcd4eec051d1122f438b94d Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 29 Jun 2026 23:06:09 -0400 Subject: [PATCH 24/31] Edited my input validation logic so that all the output matrices are empty, even if the specified number of rows and columns is not 0 --- rust/lance-linalg/src/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index c2b09dc6ebf..08ce750b681 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -153,7 +153,7 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { println!("Error: Matrix must not contain null or infinite entries."); } - if m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { + if a.len() == 0 || m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { return (vec![], vec![], vec![]); } From b67fd6ff67e68888e4500759b5d426fe7159a2c2 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 29 Jun 2026 23:49:28 -0400 Subject: [PATCH 25/31] Included comments explaining my code --- rust/lance-linalg/src/svd.rs | 84 +++++++++++++++++++++++++++++++++++- 1 file changed, 83 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 08ce750b681..3972dd91eed 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -1,6 +1,9 @@ +//Machine epsilon for f64:smallest value where 1.0 + EPS != 1.0 const EPS: f64 = 2.220_446_049_250_313e-16; +//Maximum number of Jacobi iteration sweeps before giving up const MAX_JACOBI_SWEEPS: usize = 200; +//Creates an nxn identity matrix (row-major flat vector) fn create_identity_matrix(n: usize) -> Vec { let mut m = vec![0f64; n * n]; for i in 0..n { @@ -9,17 +12,23 @@ fn create_identity_matrix(n: usize) -> Vec { m } +//Applies a Givens rotation to a symmetric matrix S in-place +//Computes S ← G^T * S * G where G is the rotation in the (p,q) plane fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { + //Reads the current values of the 2x2 submatrix at (p,p), (q,q), and (p,q) let s_pp = s[p * n + p]; let s_qq = s[q * n + q]; let s_pq = s[p * n + q]; + //Updates the 2x2 submatrix using the closed-form Jacobi rotation formulas s[p * n + p] = c * c * s_pp - 2.0 * sn * c * s_pq + sn * sn * s_qq; s[q * n + q] = sn * sn * s_pp + 2.0 * sn * c * s_pq + c * c * s_qq; s[p * n + q] = 0.0; s[q * n + p] = 0.0; + //Updates all other rows and columns that interact with p or q for r in 0..n { + //Skips the 2x2 block already handles above if r == p || r == q { continue; } @@ -28,13 +37,18 @@ fn jacobi_rotate(s: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { let new_rp = c * s_rp - sn * s_rq; let new_rq = sn * s_rp + c * s_rq; + //Updates both (r,p) and (p,r) to maintain symmetry s[r * n + p] = new_rp; s[p * n + r] = new_rp; + + //Updates both (r,q) and (q,r) to maintain symmetry s[r * n + q] = new_rq; s[q * n + r] = new_rq; } } +//Accumulates a Givens rotation into V from the right: V ← V * G +//Rotates columns p and q of V by angle (c, sn) fn apply_givens_rotation_from_right(v: &mut [f64], n: usize, p: usize, q: usize, c: f64, sn: f64) { for r in 0..n { let vp = v[r * n + p]; @@ -45,11 +59,19 @@ fn apply_givens_rotation_from_right(v: &mut [f64], n: usize, p: usize, q: usize, } } +//Jacobi eigenvalue algorithm: Decomposes a symmetric matrix A into V * diag(eigenvalues) * V^T +//by repeatedly applying Givens rotation to zero out off-diagonal entries. +//Input: a - summetric nxn matrix (row-major flat vector) +//Output: (eigenvalues, V) where the columns of V are the eigenvectors fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { + //Copy of the matrix; this copy will be diagonalized in-place. let mut s = a.to_vec(); + //Accumulates the product of all Givens rotations let mut v = create_identity_matrix(n); for _ in 0..MAX_JACOBI_SWEEPS { + //Finds the largest off-diagonal entry s[p,q] + //This is the entry we will zero out in this iteration sweep. let mut max_val = 0.0f64; let mut p = 0; let mut q = 1; @@ -63,10 +85,14 @@ fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { } } } + + //If all off-diagonal entries are negligibly small, the matrix is + //diagonal and the eigenvalues are on the diagonal. if max_val < EPS * 1e4 { break; } + //Computes the Jacobi rotation angle theta that zeros out s[p,q] let s_pq = s[p * n + q]; let diff = s[q * n + q] - s[p * n + p]; let theta = if diff.abs() < EPS { @@ -76,14 +102,22 @@ fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { }; let (sn, c) = theta.sin_cos(); + //Applies the Givens rotation: S ← G^T * S * G + //This zeroes out s[p,q] and s[q,p] while updating the rest of the matrix. jacobi_rotate(&mut s, n, p, q, c, sn); + + //Accumulates the rotation into V: V ← V * G + //At convergence, V's columns are the eigenvectors of the original matrix. apply_givens_rotation_from_right(&mut v, n, p, q, c, sn); } + //Extract eigenvalues from the diagonal of the now-diagonalized matrix let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); (eigenvalues, v) } +//Computes C = A^T A where A is a mxn row-major flat vector +//and C is an nxn row-major flat vector fn compute_ata(a: &[f64], m: usize, n: usize) -> Vec { let mut c = vec![0f64; n * n]; for i in 0..n { @@ -97,9 +131,14 @@ fn compute_ata(a: &[f64], m: usize, n: usize) -> Vec { c } +//Modified Gram-Schmidt orthonormalization +//Takes a list of column vectors and makes them orthonormal +//If a column is zero (e.g., from a zero singular value), it is +//replaced with a standard basis vector orthogonal to all prior columns. fn gram_schmidt(cols: &mut Vec>) { let ncols = cols.len(); for i in 0..ncols { + //Subtracts the projection of cols[i] onto each already-orthonormal column for j in 0..i { let dot: f64 = cols[i].iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); let cj = cols[j].clone(); @@ -108,15 +147,19 @@ fn gram_schmidt(cols: &mut Vec>) { } } + //Normalizes the resulting vector let norm = cols[i].iter().map(|&x| x * x).sum::().sqrt(); if norm > EPS { for x in cols[i].iter_mut() { *x /= norm; } } else { + //If the column is zero, finds a standard basis vector not in the span yet let dim = cols[i].len(); 'search: for k in 0..dim { let mut e = vec![0.0f64; dim]; + //Tries the k-th standard basis vector and + //orthogonalize it against all previous columns e[k] = 1.0; for j in 0..i { let dot: f64 = e.iter().zip(cols[j].iter()).map(|(&a, &b)| a * b).sum(); @@ -127,6 +170,7 @@ fn gram_schmidt(cols: &mut Vec>) { } let n2 = e.iter().map(|&x| x * x).sum::().sqrt(); if n2 > EPS { + //Found a valid replacement vector. Normalizes and uses it. for x in e.iter_mut() { *x /= n2; } @@ -138,30 +182,51 @@ fn gram_schmidt(cols: &mut Vec>) { } } +//Computes y = A * x where A is an mxn row-major flat vector, +//x is a vector of length n, and y is a vector of length m fn multiply_a_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() } +//Implements the classical eigendecomposition singular value decomposition (SVD) +//algorithm in Rust without using OpenBLAS as a dependency +//Decomposes an mxn matrix A into U * diag(sigma) * V^T +//Output: (U, sigma, V^T) where: +//U is an mxm orthogonal matrix (row-major flat vector) +//sigma is a vector of lenth min(m,n) where the singular values +//are sorted in descending order +//V^T is an nxn orthogonal matrix (row-major flat vector) +//where the rows are right singular vectors pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { + //Checks whether the input matrix A has at least 1 row and at least 1 column if a.len() == 0 || m == 0 || n == 0 { println!("Error: Matrix must have at least 1 row and at least 1 column."); } + //Checks whether the data length of the input matrix A matches the + //product of the specified number of rows and number of columns if a.len() != m * n { println!("Error: Data length of matrix must match the product of the specified number of rows and number of columns."); } + //Checks whether the input matrix A contains null or infinite entries if !a.iter().all(|x| x.is_finite()) { println!("Error: Matrix must not contain null or infinite entries."); } - + //If any of the three input validation checks above fails, + //all three output matrices are empty. if a.len() == 0 || m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { return (vec![], vec![], vec![]); } + //Step 1: Forms A^T A (nxn symmetric positive semi-definite matrix) let ata = compute_ata(a, m, n); + //Step 2: Eigendecomposes A^T A into eigenvalues λ and + //eigenvectors V using Jacobi iteration let (eigenvalues, v) = jacobi_eigen(&ata, n); + //k is the number of singular values. let k = m.min(n); + //Sorts all n eigenvalues in descending order let mut order: Vec = (0..n).collect(); for i in 1..n { let mut j = i; @@ -171,29 +236,44 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { } } + //Step 3: Computes singular values σᵢ = √λᵢ for the top k eigenvalues let mut sigma: Vec = order[..k] .iter() .map(|&i| if eigenvalues[i] > 0.0 { eigenvalues[i].sqrt() } else { 0.0 }) .collect(); + //Step 4: Computes the left singular vectors uᵢ = (1/σᵢ) * A * vᵢ + //Builds U column by column using the sorted eigenvectors let mut u_cols: Vec> = Vec::with_capacity(m); for index in 0..k { + //eᵢ is the index of the idx-th largest eigenvector in V let ei = order[index]; + //Extracts the eigenvector eᵢ (column eᵢ of V, stored as a row-major flat nxn vector) let vi: Vec = (0..n).map(|r| v[r * n + ei]).collect(); + //Computes A * vᵢ to get the unnormalized left singular vector let av = multiply_a_by_vector(a, &vi, m, n); if sigma[index] > EPS * 10.0 { + //Normalizes the left singular vector to get uᵢ u_cols.push(av.iter().map(|&x| x / sigma[index]).collect()); } else { + //The placeholder column of sigma is the zero singular value + //This column is filled by the Gram-Schmidt algorithm. sigma[index] = 0.0; u_cols.push(vec![0.0; m]); } } + //Step 5: If m > k, pad U with zero columns. + //The Gram-Schmidt algorithm will fill these columns with + //orthonormal vectors that complete the basis of R^m. for _ in k..m { u_cols.push(vec![0.0; m]); } + + //Orthonormalize all columns of U gram_schmidt(&mut u_cols); + //Pack the columns of U into an mxm flat row-major matrix let mut u = vec![0f64; m * m]; for (ci, col) in u_cols.iter().enumerate() { for (ri, &value) in col.iter().enumerate() { @@ -201,6 +281,8 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { } } + //Step 6: Build V^T (nxn row-major flat vector) + //Row i of V^T = column order[i] of V = the i-th right singular vector let mut vt = vec![0f64; n * n]; for new_row in 0..n { let old_col = order[new_row]; From 40e00acaee318c78cd4435b28e5d05c6bd7e648c Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Mon, 29 Jun 2026 23:52:03 -0400 Subject: [PATCH 26/31] Removed the tests from this file to tidy this file up --- rust/lance-linalg/src/svd.rs | 284 ----------------------------------- 1 file changed, 284 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 3972dd91eed..0215d3ba1e1 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -292,288 +292,4 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { } (u, sigma, vt) -} - -fn main() { - println!("========== UNIT TESTS =========="); - unit_test_1_multiply_a_by_vector(); - unit_test_2_compute_ata(); - unit_test_3_create_identity_matrix(); - - println!("\n========== INTEGRATION TESTS =========="); - integration_test_1_svd_2x3(); - integration_test_2_svd_4x4(); - - println!("\n========== MANUAL TESTS =========="); - manual_test_1_svd_3x3(); - println!(); - manual_test_2_svd_empty_matrix(); - println!(); - manual_test_3_svd_length_mismatch(); - println!(); - manual_test_4_svd_nan_entry(); - println!(); - manual_test_5_svd_infinite_entry(); - println!(); - manual_test_6_jacobi_rotate(); - println!(); - manual_test_7_apply_givens_rotation_from_right(); - println!(); - manual_test_8_jacobi_eigen(); - println!(); - manual_test_9_gram_schmidt(); - - println!("\nAll tests complete."); -} - -// ========================================================================= -// UNIT TESTS — automated -// ========================================================================= - -fn unit_test_1_multiply_a_by_vector() { - let a = vec![5.0, 7.0, -1.0, 4.0]; - let x = vec![-9.0, -2.0]; - let m = 2; - let n = 2; - let result = multiply_a_by_vector(&a, &x, m, n); - assert_eq!(result, vec![-59.0, 1.0], - "unit_test_1 FAILED: got {:?}, expected [-59.0, 1.0]", result); - println!("unit_test_1 PASSED"); -} - -fn unit_test_2_compute_ata() { - let a = vec![3.0, 2.0, -4.0, 7.0, -2.0, -1.0, 5.0, 4.0]; - let m = 2; - let n = 4; - let result = compute_ata(&a, m, n); - let expected = vec![ - 13.0, 8.0, -22.0, 13.0, - 8.0, 5.0, -13.0, 10.0, - -22.0, -13.0, 41.0, -8.0, - 13.0, 10.0, -8.0, 65.0, - ]; - for (i, (r, e)) in result.iter().zip(expected.iter()).enumerate() { - assert!( - (r - e).abs() < 1e-9, - "unit_test_2 FAILED at index {}: got {}, expected {}", i, r, e - ); - } - println!("unit_test_2 PASSED"); -} - -fn unit_test_3_create_identity_matrix() { - let n = 5; - let result = create_identity_matrix(n); - let expected = vec![ - 1.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 1.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 1.0, - ]; - assert_eq!(result, expected, - "unit_test_3 FAILED: got {:?}", result); - println!("unit_test_3 PASSED"); -} - -// ========================================================================= -// INTEGRATION TESTS — automated -// ========================================================================= - -fn integration_test_1_svd_2x3() { - let a = vec![3.0, 2.0, 2.0, 2.0, 3.0, -2.0]; - let m = 2; - let n = 3; - let (u, sigma, vt) = svd(&a, m, n); - - assert!(!sigma.is_empty(), "integration_test_1 FAILED: sigma is empty"); - assert!( - (sigma[0] - 5.0).abs() < 1e-3, - "integration_test_1 FAILED: sigma[0] = {}, expected 5.0", sigma[0] - ); - assert!( - (sigma[1] - 3.0).abs() < 1e-3, - "integration_test_1 FAILED: sigma[1] = {}, expected 3.0", sigma[1] - ); - - // : Vec added — this is the only change, asserts are unchanged - let expected_u: Vec = vec![0.7071, 0.7071, 0.7071, -0.7071]; - let expected_vt: Vec = vec![ - 0.7071, 0.7071, 0.0000, - 0.2357, -0.2357, 0.9428, - 0.6667, -0.6667, -0.3333, - ]; - for (i, (&r, &e)) in u.iter().zip(expected_u.iter()).enumerate() { - assert!( - (r.abs() - e.abs()).abs() < 1e-3, - "integration_test_1 FAILED U at index {}: got {}, expected {}", i, r, e - ); - } - for (i, (&r, &e)) in vt.iter().zip(expected_vt.iter()).enumerate() { - assert!( - (r.abs() - e.abs()).abs() < 1e-3, - "integration_test_1 FAILED V^T at index {}: got {}, expected {}", i, r, e - ); - } - println!("integration_test_1 PASSED: U = {:?}, sigma = {:?}, V^T = {:?}", u, sigma, vt); -} - -fn integration_test_2_svd_4x4() { - let a = vec![ - 7.0, -4.0, 5.0, 5.0, - 8.0, -2.0,-10.0, 1.0, - -1.0, -8.0, 9.0, 3.0, - 8.0, 7.0, -3.0, 4.0, - ]; - let m = 4; - let n = 4; - let (u, sigma, vt) = svd(&a, m, n); - - assert!(!sigma.is_empty(), "integration_test_2 FAILED: sigma is empty"); - - // : Vec added — asserts unchanged - let expected_sigma: Vec = vec![17.834861, 13.682368, 8.433152, 0.769723]; - for (i, (&r, &e)) in sigma.iter().zip(expected_sigma.iter()).enumerate() { - assert!( - (r - e).abs() < 1e-3, - "integration_test_2 FAILED sigma at index {}: got {}, expected {}", i, r, e - ); - } - - let expected_u: Vec = vec![ - 0.1103, 0.7669, 0.1056, -0.6233, - -0.5973, 0.3163, -0.7191, 0.1617, - 0.5994, 0.4553, -0.1537, 0.6402, - -0.5214, 0.3233, 0.6695, 0.4189, - ]; - let expected_vt: Vec = vec![ - -0.4921, -0.4312, 0.7560, -0.0187, - 0.7330, -0.3713, 0.2777, 0.4977, - 0.0588, 0.8219, 0.5131, 0.2402, - -0.4659, -0.0249, -0.2968, 0.8332, - ]; - for (i, (&r, &e)) in u.iter().zip(expected_u.iter()).enumerate() { - assert!( - (r.abs() - e.abs()).abs() < 1e-3, - "integration_test_2 FAILED U at index {}: got {}, expected {}", i, r, e - ); - } - for (i, (&r, &e)) in vt.iter().zip(expected_vt.iter()).enumerate() { - assert!( - (r.abs() - e.abs()).abs() < 1e-3, - "integration_test_2 FAILED V^T at index {}: got {}, expected {}", i, r, e - ); - } - println!("integration_test_1 PASSED: U = {:?}, sigma = {:?}, V^T = {:?}", u, sigma, vt); -} - -// ========================================================================= -// MANUAL TESTS -// ========================================================================= - -fn manual_test_1_svd_3x3() { - let a = vec![-9.0, -5.0, -2.0, 4.0, -1.0, 6.0, 9.0, -2.0, -6.0]; - let m = 3; - let n = 3; - let (u, sigma, vt) = svd(&a, m, n); - println!("manual_test_1 U: {:?}", u); - println!("manual_test_1 sigma: {:?}", sigma); - println!("manual_test_1 V^T: {:?}", vt); - println!("Expected U ≈ [-0.6955, -0.4291, -0.5763, 0.2417, 0.6157, -0.7500, 0.6767, -0.6609, -0.3245], Expected sigma ≈ [13.5037, 8.9903, 4.5633], Expected V^T ≈ [0.9861, 0.1394, -0.0903, 0.0419, 0.3172, 0.9474, -0.1607, 0.9381, -0.3070]"); -} - -fn manual_test_2_svd_empty_matrix() { - let a: Vec = vec![]; - let m = 0; - let n = 0; - let (u, sigma, vt) = svd(&a, m, n); - println!("manual_test_2 U: {:?}", u); - println!("manual_test_2 sigma: {:?}", sigma); - println!("manual_test_2 V^T: {:?}", vt); - println!("Expected: all empty vectors, error message: 'Error: Matrix must have at least 1 row and at least 1 column.'"); -} - -fn manual_test_3_svd_length_mismatch() { - let a = vec![9.0, -8.0, 3.0, -1.0]; - let m = 1; - let n = 2; - let (u, sigma, vt) = svd(&a, m, n); - println!("manual_test_3 sigma: {:?}", sigma); - println!("manual_test_3 U: {:?}", u); - println!("manual_test_3 V^T: {:?}", vt); - println!("Expected: all empty vectors, error message: 'Error: Data length of matrix must match the product of the specified number of rows and number of columns.'"); -} - -fn manual_test_4_svd_nan_entry() { - let a = vec![3.0, 1.0, f64::NAN, -4.0, -2.0, -1.0, 8.0, 3.0, 1.0]; - let m = 3; - let n = 3; - let (u, sigma, vt) = svd(&a, m, n); - println!("manual_test_4 sigma: {:?}", sigma); - println!("manual_test_4 U: {:?}", u); - println!("manual_test_4 V^T: {:?}", vt); - println!("Expected: all empty vectors, error message: 'Error: Matrix must not contain null or infinite entries.'"); -} - -fn manual_test_5_svd_infinite_entry() { - let a = vec![9.0, f64::INFINITY, -1.0, 8.0]; - let m = 2; - let n = 2; - let (u, sigma, vt) = svd(&a, m, n); - println!("manual_test_5 sigma: {:?}", sigma); - println!("manual_test_5 U: {:?}", u); - println!("manual_test_5 V^T: {:?}", vt); - println!("Expected: all empty vectors, error message: 'Error: Matrix must not contain null or infinite entries.'"); -} - -fn manual_test_6_jacobi_rotate() { - let mut a = vec![4.0, 2.0, 2.0, 3.0]; - let p = 0; - let q = 1; - let s_pq: f64 = a[p * 2 + q]; - let diff: f64 = a[q * 2 + q] - a[p * 2 + p]; - let theta: f64 = 0.5 * (2.0 * s_pq / diff).atan(); - let (sn, c) = theta.sin_cos(); - jacobi_rotate(&mut a, 2, p, q, c, sn); - println!("manual_test_6 after jacobi_rotate: {:?}", a); - println!("Off-diagonal a[0,1] = {:.6} (expected ~0.0)", a[p * 2 + q]); - println!("Trace = {:.6} (expected 7.0)", a[0] + a[3]); -} - -fn manual_test_7_apply_givens_rotation_from_right() { - let mut b = vec![1.0, 0.0, 0.0, 1.0]; - let angle: f64 = std::f64::consts::FRAC_PI_4; - let c = angle.cos(); - let s = angle.sin(); - apply_givens_rotation_from_right(&mut b, 2, 0, 1, c, s); - println!("manual_test_7 after rotation: {:?}", b); - println!("Expected: [0.7071, -0.7071, 0.7071, 0.7071]"); -} - -fn manual_test_8_jacobi_eigen() { - let a = vec![4.0, 1.0, 1.0, 3.0]; - let (eigenvalues, eigenvectors) = jacobi_eigen(&a, 2); - println!("manual_test_8 eigenvalues: {:?}", eigenvalues); - println!("manual_test_8 eigenvectors: {:?}", eigenvectors); - println!("Expected eigenvalues ≈ [4.618, 2.382] (any order)"); -} - -fn manual_test_9_gram_schmidt() { - let mut cols = vec![ - vec![1.0, 1.0, 0.0], - vec![1.0, 0.0, 1.0], - vec![0.0, 1.0, 1.0], - ]; - gram_schmidt(&mut cols); - println!("manual_test_9 after gram_schmidt:"); - for (i, col) in cols.iter().enumerate() { - let norm: f64 = col.iter().map(|&x| x * x).sum::().sqrt(); - println!(" col[{}] = {:?} norm = {:.6}", i, col, norm); - } - let dot01: f64 = cols[0].iter().zip(cols[1].iter()).map(|(&a, &b)| a * b).sum(); - let dot02: f64 = cols[0].iter().zip(cols[2].iter()).map(|(&a, &b)| a * b).sum(); - let dot12: f64 = cols[1].iter().zip(cols[2].iter()).map(|(&a, &b)| a * b).sum(); - println!("dot(col0,col1) = {:.6} (expected ~0.0)", dot01); - println!("dot(col0,col2) = {:.6} (expected ~0.0)", dot02); - println!("dot(col1,col2) = {:.6} (expected ~0.0)", dot12); } \ No newline at end of file From 49f8d5fab912ddfe70d8aa0cd1831242a728fd5c Mon Sep 17 00:00:00 2001 From: SarahNasser576 <149287747+SarahNasser576@users.noreply.github.com> Date: Fri, 17 Jul 2026 17:18:59 -0400 Subject: [PATCH 27/31] Normalized the YAML list indentation in pre-commit-config.yaml Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> --- pre-commit-config.yaml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pre-commit-config.yaml b/pre-commit-config.yaml index cbc3411268c..20cdae104a8 100644 --- a/pre-commit-config.yaml +++ b/pre-commit-config.yaml @@ -1,11 +1,11 @@ repos: -- repo: https://github.com/pre-commit/pre-commit-hooks + - repo: https://github.com/pre-commit/pre-commit-hooks rev: v2.3.0 hooks: - - id: check-yaml - - id: end-of-file-fixer - - id: trailing-whitespace -- repo: https://github.com/psf/black + - id: check-yaml + - id: end-of-file-fixer + - id: trailing-whitespace + - repo: https://github.com/psf/black rev: 22.10.0 hooks: - - id: black \ No newline at end of file + - id: black \ No newline at end of file From 2877427f39c44d1e219ff309769f159cacfee9c7 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Fri, 17 Jul 2026 17:38:41 -0400 Subject: [PATCH 28/31] Added real rustdoc comments with an example for the public API svd --- rust/lance-linalg/src/svd.rs | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 0215d3ba1e1..4bba8a288a5 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -188,15 +188,20 @@ fn multiply_a_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { (0..m).map(|i| (0..n).map(|j| a[i * n + j] * x[j]).sum()).collect() } -//Implements the classical eigendecomposition singular value decomposition (SVD) -//algorithm in Rust without using OpenBLAS as a dependency -//Decomposes an mxn matrix A into U * diag(sigma) * V^T -//Output: (U, sigma, V^T) where: -//U is an mxm orthogonal matrix (row-major flat vector) -//sigma is a vector of lenth min(m,n) where the singular values -//are sorted in descending order -//V^T is an nxn orthogonal matrix (row-major flat vector) -//where the rows are right singular vectors +/// Computes the SVD of an `m x n` row-major matrix `a`, decomposing it into +/// `U * diag(sigma) * V^T`. +/// +/// Returns `(u, sigma, vt)` where `u` is `m x m` (row-major), `sigma` has length +/// `min(m, n)` sorted descending, and `vt` is `n x n` (row-major, rows are right +/// singular vectors). +/// +/// # Example +/// ``` +/// let (u, sigma, vt) = svd(&[1.0, 0.0, 0.0, 1.0], 2, 2)?; +/// assert_eq!(u, vec![1.0, 0.0, 0.0, 1.0]); +/// assert_eq!(sigma, vec![1.0, 1.0]); +/// assert_eq!(vt, vec![1.0, 0.0, 0.0, 1.0]); +/// ``` pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { //Checks whether the input matrix A has at least 1 row and at least 1 column if a.len() == 0 || m == 0 || n == 0 { From 0f431bf860f3896ec20571a159cbfea992758ca3 Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Fri, 17 Jul 2026 18:06:22 -0400 Subject: [PATCH 29/31] Returned a Result with typed errors for invalid input to svd() function --- rust/lance-linalg/src/svd.rs | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 4bba8a288a5..5c04aee50e9 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -1,3 +1,6 @@ +//Helpful for printing error messages to determine whether +//the input (matrix, m, n) for the svd() function is valid +use lance_core::{Error, Result}; //Machine epsilon for f64:smallest value where 1.0 + EPS != 1.0 const EPS: f64 = 2.220_446_049_250_313e-16; //Maximum number of Jacobi iteration sweeps before giving up @@ -202,24 +205,26 @@ fn multiply_a_by_vector(a: &[f64], x: &[f64], m: usize, n: usize) -> Vec { /// assert_eq!(sigma, vec![1.0, 1.0]); /// assert_eq!(vt, vec![1.0, 0.0, 0.0, 1.0]); /// ``` -pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { +pub fn svd(a: &[f64], m: usize, n: usize) -> Result<(Vec, Vec, Vec)> { //Checks whether the input matrix A has at least 1 row and at least 1 column - if a.len() == 0 || m == 0 || n == 0 { - println!("Error: Matrix must have at least 1 row and at least 1 column."); + if a.is_empty() || m == 0 || n == 0 { + return Err(Error::invalid_input(format!( + "svd: matrix must have >=1 row and column (m={m}, n={n}, len={})", + a.len() + ))); } //Checks whether the data length of the input matrix A matches the //product of the specified number of rows and number of columns if a.len() != m * n { - println!("Error: Data length of matrix must match the product of the specified number of rows and number of columns."); + return Err(Error::invalid_input(format!( + "svd: data length {} != m*n ({m}*{n})", a.len() + ))); } //Checks whether the input matrix A contains null or infinite entries if !a.iter().all(|x| x.is_finite()) { - println!("Error: Matrix must not contain null or infinite entries."); - } - //If any of the three input validation checks above fails, - //all three output matrices are empty. - if a.len() == 0 || m == 0 || n == 0 || a.len() != m * n || !a.iter().all(|x| x.is_finite()) { - return (vec![], vec![], vec![]); + return Err(Error::invalid_input( + "svd: matrix contains NaN or infinite entries" + )); } //Step 1: Forms A^T A (nxn symmetric positive semi-definite matrix) @@ -296,5 +301,5 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> (Vec, Vec, Vec) { } } - (u, sigma, vt) + Ok((u, sigma, vt)) } \ No newline at end of file From d6d60e4d747c1cf129a74d1e3d072bb27f70fdff Mon Sep 17 00:00:00 2001 From: Sarah Nasser Date: Fri, 17 Jul 2026 18:25:30 -0400 Subject: [PATCH 30/31] Checked that m*n, m*m, and n*n don't overflow usize --- rust/lance-linalg/src/svd.rs | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 5c04aee50e9..6693878ac00 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -213,9 +213,30 @@ pub fn svd(a: &[f64], m: usize, n: usize) -> Result<(Vec, Vec, Vec Date: Fri, 17 Jul 2026 18:57:20 -0400 Subject: [PATCH 31/31] Made Jacobi convergence aware of small-scale matrices A and actually sweep the matrix --- rust/lance-linalg/src/svd.rs | 113 +++++++++++++++++++++++------------ 1 file changed, 74 insertions(+), 39 deletions(-) diff --git a/rust/lance-linalg/src/svd.rs b/rust/lance-linalg/src/svd.rs index 6693878ac00..56014c53d83 100644 --- a/rust/lance-linalg/src/svd.rs +++ b/rust/lance-linalg/src/svd.rs @@ -3,8 +3,9 @@ use lance_core::{Error, Result}; //Machine epsilon for f64:smallest value where 1.0 + EPS != 1.0 const EPS: f64 = 2.220_446_049_250_313e-16; -//Maximum number of Jacobi iteration sweeps before giving up -const MAX_JACOBI_SWEEPS: usize = 200; +//Maximum number of Jacobi sweeps before giving up, when not +//scaled by dimension +const MAX_JACOBI_SWEEPS: usize = 100; //Creates an nxn identity matrix (row-major flat vector) fn create_identity_matrix(n: usize) -> Vec { @@ -66,57 +67,91 @@ fn apply_givens_rotation_from_right(v: &mut [f64], n: usize, p: usize, q: usize, //by repeatedly applying Givens rotation to zero out off-diagonal entries. //Input: a - summetric nxn matrix (row-major flat vector) //Output: (eigenvalues, V) where the columns of V are the eigenvectors -fn jacobi_eigen(a: &[f64], n: usize) -> (Vec, Vec) { +fn jacobi_eigen(a: &[f64], n: usize) -> Result<(Vec, Vec)> { //Copy of the matrix; this copy will be diagonalized in-place. let mut s = a.to_vec(); //Accumulates the product of all Givens rotations let mut v = create_identity_matrix(n); - for _ in 0..MAX_JACOBI_SWEEPS { - //Finds the largest off-diagonal entry s[p,q] - //This is the entry we will zero out in this iteration sweep. - let mut max_val = 0.0f64; - let mut p = 0; - let mut q = 1; - for i in 0..n { - for j in i + 1..n { - let val = s[i * n + j].abs(); - if val > max_val { - max_val = val; - p = i; - q = j; + //There is nothing to do for 1x1 matrices due to there being no off-diagonal entries. + if n <= 1 { + let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); + return Ok((eigenvalues, v)); + } + + //Uses the Frobenius norm of the input matrix as a scale reference for convergence + let norm: f64 = s.iter().map(|&x| x * x).sum::().sqrt(); + if norm < EPS { + let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); + return Ok((eigenvalues, v)); + } + + //Relative tolerance: off-diagonal entries (and the overall off-diagonal + //Frobenius norm) smaller than this, relative to the matrix's own scale, + //are converged + let relative_tolerance = EPS * norm * (n as f64); + + //Scales the sweep budget with the matrix dimension + //Larger matrices need more sweeps and individual rotations to + //fully diagonalize. + let max_sweeps = MAX_JACOBI_SWEEPS.max(n * n); + + let mut converged = false; + for _sweep in 0..max_sweeps { + //This is one full cyclic sweep. + //It visits each (p, q) pair with p < q exactly once, + //applying a rotation when the entry is not already negligible. + for p in 0..n { + for q in p + 1..n { + let s_pq = s[p * n + q]; + if s_pq.abs() <= relative_tolerance { + continue; } + + //Computes the Jacobi rotation angle theta that zeros out s[p,q] + let diff = s[q * n + q] - s[p * n + p]; + let theta = if diff.abs() < EPS { + std::f64::consts::FRAC_PI_4 + } else { + 0.5 * (2.0 * s_pq / diff).atan() + }; + let (sn, c) = theta.sin_cos(); + + //Applies the Givens rotation: S ← G^T * S * G + //This zeroes out s[p,q] and s[q,p] while updating the rest of the matrix. + jacobi_rotate(&mut s, n, p, q, c, sn); + + //Accumulates the rotation into V: V ← V * G + //At convergence, V's columns are the eigenvectors of the original matrix. + apply_givens_rotation_from_right(&mut v, n, p, q, c, sn); } } - - //If all off-diagonal entries are negligibly small, the matrix is - //diagonal and the eigenvalues are on the diagonal. - if max_val < EPS * 1e4 { + + //Recomputes the off-diagonal Frobenius norm after the full sweep to + //determine whether another sweep is needed + let mut off_diag_sq = 0.0f64; + for p in 0..n { + for q in p + 1..n { + off_diag_sq += s[p * n + q] * s[p * n + q]; + } + } + if off_diag_sq.sqrt() < relative_tolerance { + converged = true; break; } + } - //Computes the Jacobi rotation angle theta that zeros out s[p,q] - let s_pq = s[p * n + q]; - let diff = s[q * n + q] - s[p * n + p]; - let theta = if diff.abs() < EPS { - std::f64::consts::FRAC_PI_4 - } else { - 0.5 * (2.0 * s_pq / diff).atan() - }; - let (sn, c) = theta.sin_cos(); - - //Applies the Givens rotation: S ← G^T * S * G - //This zeroes out s[p,q] and s[q,p] while updating the rest of the matrix. - jacobi_rotate(&mut s, n, p, q, c, sn); - - //Accumulates the rotation into V: V ← V * G - //At convergence, V's columns are the eigenvectors of the original matrix. - apply_givens_rotation_from_right(&mut v, n, p, q, c, sn); + //Checks whether the Jacobi eigensolver converged within the + //maximum number of sweeps + if !converged { + return Err(Error::invalid_input(format!( + "svd: Jacobi eigensolver failed to converge within {max_sweeps} sweeps for a {n}x{n} matrix" + ))); } - //Extract eigenvalues from the diagonal of the now-diagonalized matrix + //Extracts eigenvalues from the diagonal of the now-diagonalized matrix let eigenvalues: Vec = (0..n).map(|i| s[i * n + i]).collect(); - (eigenvalues, v) + Ok((eigenvalues, v)) } //Computes C = A^T A where A is a mxn row-major flat vector