2a:null 2b:["$","title",null,{"children":"Least Squares Curve Fit with the QR Decomposition | Ricky Reusser | Observable"}] 30:T949,function fitLeastSquares(x, y, basisFns) { if (x.length !== y.length) { throw new Error("Incompatible input vectors, x.length != y.length"); } if (x.length < basisFns.length) { throw new Error("Insufficient data, x.length < basisFns.length"); } const numBasis = basisFns.length; const N = x.length; // Design matrix A (row-major flat array of size N x M) const A = Array(N * numBasis).fill(0); for (let j = 0; j < numBasis; j++) { for (let i = 0; i < N; i++) { A[i * numBasis + j] = basisFns[j](x[i]); } } function column(A, rows, cols, j) { const col = Array(rows).fill(0); for (let i = 0; i < rows; i++) col[i] = A[i * cols + j]; return col; } function dot(a, b) { let sum = 0; for (let i = 0; i < a.length; i++) sum += a[i] * b[i]; return sum; } // QR Decomposition (Gram-Schmidt) function gramSchmidtQRDecomp(A, N, M) { const Q = Array(N * M).fill(0); const R = Array(M * M).fill(0); for (let j = 0; j < M; j++) { const v = column(A, N, M, j); for (let i = 0; i < j; i++) { const qi = column(Q, N, M, i); const rij = dot(qi, v); R[i * M + j] = rij; for (let k = 0; k < N; k++) { v[k] -= rij * qi[k]; } } const rjj = Math.sqrt(dot(v, v)); R[j * (M + 1)] = rjj; for (let k = 0; k < N; k++) { Q[k * M + j] = v[k] / rjj; } } return { Q, R }; } // Back substitution function backSubstitute(R, Qt_y, M) { const beta = Array(M).fill(0); for (let i = M - 1; i >= 0; i--) { let sum = Qt_y[i]; for (let j = i + 1; j < M; j++) { sum -= R[i * M + j] * beta[j]; } beta[i] = sum / R[i * (M + 1)]; } return beta; } const { Q, R } = gramSchmidtQRDecomp(A, N, numBasis); // Compute Q^T y const Qt_y = Array(numBasis).fill(0); for (let i = 0; i < N; i++) { for (let j = 0; j < numBasis; j++) { Qt_y[j] += Q[i * numBasis + j] * y[i]; } } const coefficients = backSubstitute(R, Qt_y, numBasis); let sumSquares = 0.0; for (let i = 0; i < N; i++) { let Ax_i = 0; for (let j = 0; j < numBasis; j++) { Ax_i += A[i * numBasis + j] * coefficients[j]; } const residual = y[i] - Ax_i; sumSquares += residual ** 2; } return { coefficients, error: Math.sqrt(sumSquares / N) }; }31:T609,signature(fitLeastSquares, { description: md` Use a [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition) to perform a least squares fit of arbitrary basis functions to a set of ${tex`(x_i, y_i)`} input data. **Import**: \`\`\`js import { fitLeastSquares } from '${await getPinnedSlug( "@rreusser/least-squares-curve-fit" )}' \`\`\` **Parameters**: - \`x: number[]\`: ${tex`x`}-coordinate of input data - \`y: number[]\`: ${tex`y`}-coordinate of input data - \`basisFns: Array<(number) => number>\`: list of basis functions **Returns**: - \`coefficients: number[]\`: list of coefficients for the respective basis functions. - \`error: number\`: root mean square error (RMSE) **Notes**: - The length of arrays \`x\` and \`y\` must be equal. - Length of input data must be greater than or equal to the number of basis functions.`, example: [ ` // Fit a quadratic polynomial to y = 1 + x² fitLeastSquares( [0, 1, 2, 3, 4], [1, 2, 5, 10, 17], [(x) => 1, (x) => x, (x) => x**2] ) // => { coefficients: [ 1, 0, 1 ], error: 0 }` ], tests: { "basic curve fit": (assert) => { const fit = fitLeastSquares( [0, 1, 2, 3, 4], [1, 2, 5, 10, 17], [(x) => 1, (x) => x, (x) => x ** 2] ); const expected = [1, 0, 1]; assert(Math.abs(fit.coefficients[0] - expected[0]) < 1e-8); assert(Math.abs(fit.coefficients[1] - expected[1]) < 1e-8); assert(Math.abs(fit.coefficients[2] - expected[2]) < 1e-8); assert(Math.abs(fit.error - 0 < 1e-8)); } } })32:[] 2e:{"id":"19749e186dfba23c","slug":"least-squares-curve-fit","title":"Least Squares Curve Fit with the QR Decomposition","thumbnail":"2f6b21d1b419e29f4900708a6106b44846f8b3a9a61bb6db0aeb20309597bd23","access_level":"public","create_time":"$D2025-07-20T03:49:01.161Z","update_time":"$D2025-07-21T22:35:13.006Z","view_time":null,"first_public_time":"$D2025-07-20T04:16:39.451Z","has_unpublished_edits":null,"fork":false,"stars":3,"starred":false,"owner":{"type":"team","id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/49c18a2ac0092176d4ea8348a1dd538e0b460b086b89d602b8d21ffe0821d92e","tier":"starter_2024"},"version":355,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null}],"members":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null,"role":"editor"},{"type":"team","id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/49c18a2ac0092176d4ea8348a1dd538e0b460b086b89d602b8d21ffe0821d92e","tier":"starter_2024","role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":{"id":"dcd095ad05c2b722","slug":"writeups","title":"Fun Notebooks","owner":{"id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser"}},"type":"notebook","body":{"titling":"auto","title":"Least Squares Curve Fit with the QR Decomposition","cells":[{"id":0,"value":"# Least Squares Curve Fit with the QR Decomposition","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"","database":"$undefined"},{"id":4,"value":"$30","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":219,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":281,"value":"## Example: Fit a fifth-order polynomial to ${tex`\\sin(x)`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"","database":"$undefined"},{"id":302,"value":"N = 51","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":45,"value":"x = Array.from({ length: N }, (_, i) => (Math.PI * 2 * i) / (N - 1))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":87,"value":"y = x.map(x => Math.sin(x))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":78,"value":"basisFns = [\n (x) => x ** 0,\n (x) => x ** 1,\n (x) => x ** 2,\n (x) => x ** 3,\n (x) => x ** 4,\n (x) => Math.exp(x) // Doesn't have to be a polynomial!\n]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":60,"value":"fit = fitLeastSquares(x, y, basisFns)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":82,"value":"{\n const N = 101;\n const yFit = x.map((x) =>\n basisFns\n .slice(0, fit.coefficients.length)\n .reduce((sum, fn, i) => sum + fit.coefficients[i] * fn(x), 0)\n );\n const plotData = x.map((x, i) => ({\n x: x,\n y: y[i],\n yFit: yFit[i]\n }));\n\n return Plot.plot({\n grid: true,\n marks: [\n Plot.dot(plotData, {\n x: \"x\",\n y: \"y\",\n r: 2.5,\n fill: \"black\",\n stroke: \"none\"\n }),\n Plot.line(plotData, {\n x: \"x\",\n y: \"yFit\",\n stroke: \"red\",\n strokeWidth: 2\n })\n ],\n marginLeft: 40,\n marginBottom: 40\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"},{"id":215,"value":"import { signature, getPinnedSlug } from \"@mootari/signature@553\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W32","stdlib":"1"}} 2f:"$2e:body"