30:T671,md`# Linear Algebra via Eigen and WebAssembly: A proof of concept This notebook is heavily inspired by Rohan Sawhney's [linear-algebra-js](https://rohan-sawhney.github.io/linear-algebra-js/) and explores what it means to use the C++ [Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) library ergonomically in Observable. It's compiled to JavaScript using the [\`embind\`](https://emscripten.org/docs/porting/connecting_cpp_and_javascript/embind.html) tool of [Emscripten](https://emscripten.org/index.html). ***This notebook represents an experiment which is certainly not complete and is subject to changes and extensions or just throwing it out and trying a different strategy. You certainly should not rely on it, or at the very least not without pinning the version.*** I'd love to either tighten this up and publish it thoughtfully or just throw it all out and write some more Observable-focused glue code around the [linear-algebra-js](https://rohan-sawhney.github.io/linear-algebra-js/) library. Your input is welcome! Would love to discuss on the Observable forum. To start though, here is a wish list of features I'd still like to include: - Vector data types (as opposed to row or column matrices) - Sparse matrices - A double precision alternative - Clean up the semantics of memory management - Figure out how to distribute it cleanly via npm or github - Rows, columns, blocks, and diagonals are not strictly matrices but are used to *modify* their matrices. It'd be nice if it were clean and coherent how to *use* (or at least cast) them as matrices. I'm tempted to throw them out entirely as unnecessarily complicated.`31:T422,md`## Contents - [Eigen via \`embind\`](#eigenViaEmbind) - [Basic usage](#basicUsage) - [Avoiding leaks](#avoidingLeaks) - [Direct inversion](#directInversion) - [Matrices from JavaScript data](#matricesFromJS) - [Decompositions](#decompositions) - [Partial-pivot LU](#partialPivotLUDecomposition) - [Full-pivot LU](#fullPivotLUDecomposition) - [Cholesky LLT](#choleskyDecomposition) - [Robust Cholesky LDLT](#robustCholeskyDecomposition) - [Householder QR](#householderQRDecomposition) - [Column-pivot Householder QR](#columnPivotQRDecomposition) - [Full-pivot Householder QR](#fullPivotQRDecomposition) - [Bidiagonal Divide and Conquer SVD](#BDCSVDDecomposition) - [Self Adjoint Eigensolver](#SelfAdjointEigenSolver) - [Interfacing with NDarrays](#interfacingWithNdarray) - [Displaying matrices](#displaying) - [Block and diagonal methods](#blockAndDiag) - [Diagonals](#diagonals) - [Blocks](#blocks) - [Rows and columns](#rowsCols) - [Norms](#norms) - [Advanced memory management](#advancedMemoryManagement) - [License](#LICENSE) `32:T764,function manageInvalidation(Library, wrappedMethods) { var globalInvalidation = null; Library.ensure = function ensure(invalidation, callback) { var error; try { let resolveTmp; const tmp = new Promise((resolve, reject) => { resolveTmp = resolve; }); globalInvalidation = invalidation; var returnValue = callback.bind(this)(tmp); resolveTmp(); return returnValue; } catch (e) { error = e; } finally { globalInvalidation = null; } if (error) throw error; }; function ensureInvalidation(obj) { return function(name) { var fn = obj[name]; if (!fn) throw new Error(name + ' is not defined'); var nArgs = fn.length; var fnName = fn .toString() .match(/function ([^\(]*)/)[1] .split('$'); var wrappedFn = function() { let ret; for (var i = 0, args = []; i < nArgs; i++) args[i] = arguments[i]; ret = fn.apply(this, args); var invalidation = arguments[nArgs]; var i = invalidation || globalInvalidation; if (!i || i.then === undefined) { throw new Error( `To safely manage memory, you must either wrap this operation inside an \`.ensure(invalidation, () => ...)\` block or pass \`invalidation\` as argument ${nArgs + 1} of ${fnName.join('.')}.` ); } i.then(() => { if (!ret.isDeleted()) ret.delete(); }); return ret; }; Object.defineProperty(wrappedFn, "name", { value: fnName.join('::') }); obj[name] = wrappedFn; }; } Object.keys(wrappedMethods).forEach(c => { (wrappedMethods[c].class || []).forEach(ensureInvalidation(Library[c])); (wrappedMethods[c] || []).prototype.forEach( ensureInvalidation(Library[c].prototype) ); }); return Library; }33:T6b9,managedMethods = ({ Matrixf: { prototype: [ 'inverse', 'transpose', 'diagonal', 'block', 'map', 'topLeftCorner', 'topRightCorner', 'bottomLeftCorner', 'bottomRightCorner', 'topRows', 'bottomRows', 'leftCols', 'rightCols', 'mul', 'add', 'sub', 'muls', 'adds', 'subs', 'ldlt', 'llt', 'svd', 'lu', 'qr', 'fullPivLU', 'partialPivLU', 'householderQR', 'colPivHouseholderQR', 'fullPivHouseholderQR', 'bdcSVD', 'selfAdjointEigenSolver' ], class: [ 'fromArray', 'fromDiagonal', 'fromBlock', 'fromRowMajorBlock', 'fromRow', 'fromCol', 'columnVector', 'rowVector', 'from', 'ones', 'zeros', 'random', 'identity' ] }, FullPivLUf: { prototype: ['matrixU', 'matrixL', 'permutationQ', 'permutationP', 'solve'] }, PartialPivLUf: { prototype: ['matrixU', 'matrixL', 'permutationP', 'solve'] }, LLTf: { prototype: ['matrixL', 'solve'] }, LDLTf: { prototype: ['matrixU', 'matrixL', 'transpositionsP', 'vectorD', 'solve'] }, HouseholderQRf: { prototype: ['householderQ', 'solve'] }, ColPivHouseholderQRf: { prototype: [ 'colsPermutation', 'householderQ', 'solve', 'matrixR', 'matrixQR', 'inverse' ] }, FullPivHouseholderQRf: { prototype: ['colsPermutation', 'matrixQ', 'solve', 'inverse', 'matrixQR'] }, BDCSVDf: { prototype: ['solve', 'matrixU', 'matrixV', 'singularValues'] }, SelfAdjointEigenSolverf: { prototype: ['eigenvalues', 'eigenvectors'] } })34:T432,LICENSE = "mit" /* Copyright 2019 Ricky Reusser Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */35:[{"id":"c9e03650eab306fafcefa082e91d4c35c900b61b1f80bc4c8bb8340fe32d2c7069f8268cd3647d9e79f394a4531eabdeb48e3716454b87b3d33361098c3c3bf0","href":"https://static.observableusercontent.com/files/c9e03650eab306fafcefa082e91d4c35c900b61b1f80bc4c8bb8340fe32d2c7069f8268cd3647d9e79f394a4531eabdeb48e3716454b87b3d33361098c3c3bf0","name":"tweet.png","type":"image/png","lastModified":1576012874273.631,"size":140599},{"id":"ba77744394489c0c77a0fad7e63139152dc9342bf97b0d7e6523647f58072129ede43a8b98aadebf389b82e7a9043a6e68675f42e8e405b3d8ea7c570786fe96","href":"https://static.observableusercontent.com/files/ba77744394489c0c77a0fad7e63139152dc9342bf97b0d7e6523647f58072129ede43a8b98aadebf389b82e7a9043a6e68675f42e8e405b3d8ea7c570786fe96","name":"library.min@45.js","type":"application/javascript","lastModified":1575928088471.02,"size":550343},{"id":"bf98a8e54bbd68d8077aabc19c76d9e92692911c7a766610fb906b553c28ed38441d839389807dda932570ca629aa6844746fe25baa08105a24b7f5f91e67b6a","href":"https://static.observableusercontent.com/files/bf98a8e54bbd68d8077aabc19c76d9e92692911c7a766610fb906b553c28ed38441d839389807dda932570ca629aa6844746fe25baa08105a24b7f5f91e67b6a","name":"library@7.cpp","type":"application/octet-stream","lastModified":1575928096062.329,"size":28162},{"id":"737d01511a3e0edd9811ed551ebbf5a11ee87d2e6c2acff643adc2df1d03d2a01c5d1afcd42ad6bd93f9e5f6627a98f5db6c98282a72407c0ab576f2c99657a8","href":"https://static.observableusercontent.com/files/737d01511a3e0edd9811ed551ebbf5a11ee87d2e6c2acff643adc2df1d03d2a01c5d1afcd42ad6bd93f9e5f6627a98f5db6c98282a72407c0ab576f2c99657a8","name":"eigen.cpp","type":"application/octet-stream","lastModified":1589665426012.545,"size":29023},{"id":"5b20fe80b95f6ce95e7d6ca23338bd46d2c656deacb7a0b61ef9c945644a71fba808d7aa37cac1c0f0ccd089195c8f345727dc62c77d72e257f5f85be6ba373a","href":"https://static.observableusercontent.com/files/5b20fe80b95f6ce95e7d6ca23338bd46d2c656deacb7a0b61ef9c945644a71fba808d7aa37cac1c0f0ccd089195c8f345727dc62c77d72e257f5f85be6ba373a","name":"eigen@3.js","type":"application/javascript","lastModified":1589667267857.77,"size":51354},{"id":"d762068f83d5f492c3bd2ee843e21f0bcc31a8d31b3fffa7622fee84b5f1feda6d69a1f71f42109d960fe29ff86446bb2b8e7d321c15b1d7b65415c5528c8559","href":"https://static.observableusercontent.com/files/d762068f83d5f492c3bd2ee843e21f0bcc31a8d31b3fffa7622fee84b5f1feda6d69a1f71f42109d960fe29ff86446bb2b8e7d321c15b1d7b65415c5528c8559","name":"eigen@3.wasm","type":"application/wasm","lastModified":1589667270946.013,"size":388026},{"id":"57c8c5c0ed42874ee5175a92c21fb707f09e6eb5297f62dd7bb9c21bffe47f442e91ce17138406715813a6e88c36aae00636155dd90a22795ddcc4c477d4adc8","href":"https://static.observableusercontent.com/files/57c8c5c0ed42874ee5175a92c21fb707f09e6eb5297f62dd7bb9c21bffe47f442e91ce17138406715813a6e88c36aae00636155dd90a22795ddcc4c477d4adc8","name":"eigen@1.cpp","type":"application/octet-stream","lastModified":1589680784319.452,"size":29288}] 2e:{"title":"Linear Algebra via Eigen and WebAssembly: A proof of concept","cells":[{"id":283,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":2714,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":650,"value":"viewof eigenViaEmbind = md`## Eigen via Embind\n\nThe C++ code below implements bindings with which we can utilize the Eigen linear algebra library in JavaScript. We attach the compiled result to this notebook as a file attachment, then load it and perform some additional wrapping in JavaScript to obtain the API we desire. Overall the required wrapping is very minimal and is currently limited to two features:\n\n1. enforcing memory management of allocated results via \\`invalidation\\`\n2. augmenting matrices with a \\`show()\\` method`","pinned":false,"hidden":false,"mode":"ojs"},{"id":222,"value":"{\n var cpp = await FileAttachment(\"eigen@1.cpp\").text();\n return md`\n

eigen.cpp

\n${md`\n\\`\\`\\`cpp\n${cpp}\n\\`\\`\\`\n`}\n

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2242,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":2339,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2204,"value":"Eigen = {\n window.EigenWASMPath = await FileAttachment(\"eigen@3.wasm\").url();\n var url = await FileAttachment(\"eigen@3.js\").url();\n var Initialize = await require(url);\n return new Promise((resolve, reject) => {\n Initialize().then(Eigen => {\n delete Eigen.then;\n try {\n Eigen = manageInvalidation(Eigen, managedMethods);\n\n Eigen.Matrixf.prototype.show = function() {\n for (var i = 0, args = [this]; i < arguments.length; i++)\n args[i + 1] = arguments[i];\n return showMatrix.apply(this, args);\n };\n\n resolve(Eigen);\n } catch (e) {\n reject(e);\n }\n });\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2178,"value":"Matrixf = Eigen.Matrixf","pinned":true,"hidden":false,"mode":"ojs"},{"id":1068,"value":"viewof basicUsage = md`## Some basic matrices`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1071,"value":"Matrixf.zeros([2, 2], invalidation).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":1079,"value":"Matrixf.identity([2, 2], invalidation).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":1572,"value":"Matrixf.random([2, 4], invalidation).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":550,"value":"viewof avoidingLeaks = md`## Avoiding memory leaks\n\nThe Emscripten/WASM heap and JavaScript do not play well together. Unless we explicitly deallocate matrices, we'll leak memory—potentialy a *lot*.\n\nIn this notebook, all of the functions which allocate new matrices *must* receive the built-in \\`invalidation\\` promise so that they are deleted when they go out of scope.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":558,"value":"Matrixf.identity([3, 3]).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":721,"value":"md`As long as we clean up after ourselves, we can let Observable reevaluate cells all day and not leak memory without bound.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":578,"value":"viewof matricesFromJS = md`## Creating matrices from JavaScript data\n\nJavaScript data lives in JavaScript. C++ data lives on the WebAssembly heap. And there's really no way around that. What we *can* do though is use fast copying so that it's at least as fast as it can be to transfer data across that boundary.\n\nThe code below creates native JavaScript data and then uses it to construct a large ${tex`${m} \\times ${m}`} matrix.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":672,"value":"viewof m = slider({\n min: 1,\n max: 2048,\n step: 1,\n value: 1000\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":584,"value":"data = new Float32Array(new Array(m * m).fill(0).map((d, i) => i))","pinned":true,"hidden":false,"mode":"ojs"},{"id":590,"value":"// BAD! This leaks memory since the promise never resolves, triggering cleanup!!\nMatrixf.fromArray(data, [m, m], new Promise(() => {})).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":715,"value":"{\n m;\n return md`Observe that the Emscripten heap size (\\`Eigen.HEAP8.length\\`) starts out at 16777216 bytes, but if you increase the size of the beyond a certain threshold, it will get doubled.\n\n${html`Heap size: ${Eigen.HEAP8.length / 1024 / 1024} MB`}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":531,"value":"viewof directInversion = md`## Solving ${tex`Ax = b`}\n\n${tex`Ax = b`} is *the* classic linear algebra problem whose solution is ${tex`x = A^{-1} b`}. You don't need to invert ${tex`A`} to solve for ${tex`x`}, but we'll keep it simple and do just that to start. We first construct a (symmetric) ${`A`} matrix using \\`MatrixXf.from\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2253,"value":"viewof n = slider({ min: 1, max: 10, step: 1, value: 6 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":488,"value":"viewof A = Matrixf.from(\n (i, j) => [-1, 2, -1][j - i + 1] || 0,\n [n, n],\n invalidation\n).show('A')","pinned":true,"hidden":false,"mode":"ojs"},{"id":504,"value":"viewof b = Matrixf.random([n, 1], invalidation).show('b')","pinned":true,"hidden":false,"mode":"ojs"},{"id":509,"value":"viewof x = A.inverse(invalidation)\n .mul(b, invalidation)\n .show('x = A^{-1} b')","pinned":true,"hidden":false,"mode":"ojs"},{"id":2501,"value":"md`To avoid needing to type \\`invalidation\\` more than once per cell, we create an \\`ensure\\` helper which ensures all memory allocated in the block is freed when the invalidation promise resolves. [See below](#advancedMemoryManagement) for additional details.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2234,"value":"Eigen.ensure(invalidation, () =>\n A.inverse()\n .mul(b)\n .transpose()\n).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":3122,"value":"As for the operator overloading, we can hope.\n\n${html`

`}`}`","pinned":false,"hidden":false,"mode":"md"},{"id":2594,"value":"md`## Decompositions\n\nThe Eigen documentation includes [a full list and comparison of the decompositions it implements](https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html). The section below implements them and uses them to solve ${`Ax = b`}, though in general there is a lot more you may do with them than solve a single equation.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1952,"value":"viewof partialPivotLUDecomposition = md`### Partial-pivot LU Decomposition`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2603,"value":"html`

Requirements:: Invertible
Speed:: Fast
Reliability:: Depends on condition number
Rank-revealing:: -
Also computes:: -
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1955,"value":"LU = A.partialPivLU(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1963,"value":"Eigen.ensure(invalidation, () =>\n LU.permutationP()\n .inverse()\n .mul(LU.matrixL())\n .mul(LU.matrixU())\n).show('A = P^{-1}LU')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1958,"value":"LU.solve(b, invalidation).show('b')","pinned":true,"hidden":false,"mode":"ojs"},{"id":1363,"value":"viewof fullPivotLUDecomposition = md`### Full-pivot LU Decomposition`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2614,"value":"html`

Requirements:: -
Speed:: Slow
Reliability:: Proven
Rank-revealing:: Yes
Also computes:: -
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1366,"value":"fullPivLU = A.fullPivLU(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1811,"value":"Eigen.ensure(invalidation, () =>\n fullPivLU\n .permutationP()\n .inverse()\n .mul(fullPivLU.matrixL())\n .mul(fullPivLU.matrixU())\n .mul(fullPivLU.permutationQ().inverse())\n).show('A = P^{-1}LUQ^{-1}')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1575,"value":"fullPivLU.solve(b, invalidation).show('b')","pinned":true,"hidden":false,"mode":"ojs"},{"id":1666,"value":"viewof choleskyDecomposition = md`### Cholesky Decomposition`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2682,"value":"html`\n

Requirements:: Positive definite
Speed:: Very fast
Reliability:: Depends on condition number
Rank-revealing:: -
Also computes:: -
Documentation: 🔗Link

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1772,"value":"LLT = A.llt(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1821,"value":"Eigen.ensure(invalidation, () =>\n LLT.matrixL().mul(LLT.matrixL().transpose())\n).show('A = LL^T')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1775,"value":"LLT.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":1766,"value":"viewof robustCholeskyDecomposition = md`### Robust Cholesky Decomposition`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2688,"value":"html`\n

Requirements:: Positive or negative semidefinite
Speed:: Very fast
Reliability:: Good
Rank-revealing:: -
Also computes:: -
Documentation: 🔗Link

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1669,"value":"LDLT = A.ldlt(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1913,"value":"Eigen.ensure(invalidation, () =>\n LDLT.transpositionsP()\n .mul(LDLT.matrixL())\n .mul(Matrixf.fromDiagonal(LDLT.vectorD()))\n .mul(LDLT.matrixU())\n).show('A = PLDL^{T}')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1687,"value":"LDLT.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":1987,"value":"viewof householderQRDecomposition = md`### Householder QR`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2627,"value":"html`

Requirements:: -
Speed:: Fast
Reliability:: Depends on condition number
Rank-revealing:: -
Also computes:: Orthogonalization
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1990,"value":"tex`A = QR`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1993,"value":"QR = A.householderQR(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1996,"value":"QR.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":2010,"value":"viewof columnPivotQRDecomposition = md`### Column-Pivot Householder QR`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2634,"value":"html`

Requirements:: -
Speed:: Fast
Reliability:: Good
Rank-revealing:: Yes
Also computes:: Orthogonalization
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2013,"value":"colPivQR = A.colPivHouseholderQR(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2016,"value":"colPivQR.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":2022,"value":"viewof fullPivotQRDecomposition = md`### Full-Pivot Householder QR`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2638,"value":"html`

Requirements:: -
Speed:: Slow
Reliability:: Proven
Rank-revealing:: Yes
Also computes:: Orthogonalization
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2020,"value":"fullPivQR = A.fullPivHouseholderQR(invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2025,"value":"fullPivQR.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":2040,"value":"viewof BDCSVDDecomposition = md`### Bidiagonal Divide and Conquer SVD`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2641,"value":"html`

Requirements:: -
Speed:: One of the fastest SVD algorithms
Reliability:: Excellent
Rank-revealing:: Yes
Also computes:: Singular values/vectors, least squares
Documentation: 🔗Link

\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2048,"value":"svd = A.bdcSVD(Eigen.ComputeThinU | Eigen.ComputeThinV, invalidation)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2119,"value":"svd.singularValues(invalidation).show('diag(\\\\Sigma)')","pinned":false,"hidden":false,"mode":"ojs"},{"id":2058,"value":"Eigen.ensure(invalidation, () =>\n svd\n .matrixU()\n .mul(Matrixf.fromDiagonal(svd.singularValues()))\n .mul(svd.matrixV().transpose())\n .map((i, j, value) => Math.round(value * 1000) / 1000)\n).show('A = U \\\\Sigma V^T')","pinned":false,"hidden":false,"mode":"ojs"},{"id":2096,"value":"svd.solve(b, invalidation).show('x')","pinned":true,"hidden":false,"mode":"ojs"},{"id":3197,"value":"viewof SelfAdjointEigenSolver = md`### Eigenvalues and eigenvectors`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3202,"value":"eig = Eigen.ensure(invalidation, () =>\n A.selfAdjointEigenSolver(Eigen.ComputeEigenvectors)\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":3213,"value":"viewof eigenvalues = eig.eigenvalues(invalidation).show('diga(\\\\Lambda)')","pinned":true,"hidden":false,"mode":"ojs"},{"id":3216,"value":"viewof eigenvectors = eig.eigenvectors(invalidation).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":3232,"value":"Eigen.ensure(invalidation, () =>\n eigenvectors\n .mul(Matrixf.fromDiagonal(eigenvalues))\n .mul(eigenvectors.transpose())\n .map((i, j, value) => Math.round(value * 1000) / 1000)\n).show('A = X \\\\Lambda X^T')","pinned":false,"hidden":false,"mode":"ojs"},{"id":619,"value":"viewof interfacingWithNdarray = md`## Interfacing with ndarrays\n\nIt's very easy to feed typed array views of the WebAssembly heap into either [@stdlib/ndarray](https://github.com/stdlib-js/stdlib/tree/develop/lib/node_modules/%40stdlib/ndarray/ctor) or [@scijs/ndarray](https://github.com/scijs/ndarray).\n\n***You must be careful not to retain views of the Emscripten heap since they become invalid if objects are deleted or even just if Emscripten expands the heap. The remedy is simple: you must instead acquire and use them synchronously, when you need them.***`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2981,"value":"({\n buffer: A.buffer.slice(),\n shape: A.shape,\n stride: A.stride,\n offset: A.offset\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":2563,"value":"md`### Using [\\`@scijs/ndarray\\`](https://github.com/scijs/ndarray)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"import { ndarray } from 'f181985b512ab7f4'","pinned":false,"hidden":false,"mode":"ojs"},{"id":622,"value":"ndarray(A.buffer, A.shape, A.offset, A.stride)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2567,"value":"md`### Using [\\`@stdlib/ndarray\\`](https://github.com/stdlib-js/stdlib/tree/develop/lib/node_modules/%40stdlib/ndarray/array)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":732,"value":"stdlib = require('https://unpkg.com/@stdlib/stdlib@0.0.90/dist/stdlib-flat.min.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":911,"value":"stdlib.array(A.buffer, {\n shape: A.shape,\n stride: A.stride,\n offset: A.offset,\n order: 'column-major'\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":628,"value":"md`If you need to retain data in JavaScript or do not want to modify values on the Emscripten heap, you may instead simply copy the buffer out of the WASM heap with \\`.slice()\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":634,"value":"ndarray(A.buffer.slice(), A.shape, A.stride, A.offset)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2532,"value":"viewof displaying = md`## Displaying matrices`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2768,"value":"md`You may output matrix to a string using Eigen's built-in conversion with \\`.toString()\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2772,"value":"A.toString()","pinned":true,"hidden":false,"mode":"ojs"},{"id":2774,"value":"md`For output to a nicely formatted matrix, this notebook also augments the \\`Matrixf\\` prototype with a \\`show()\\` function.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3,"value":"import { createNdarrayShow } from 'b616fe201223136a'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2300,"value":"showMatrix = {\n const showNdarray = createNdarrayShow({\n maxRows: 10,\n maxCols: 10,\n precision: 4\n });\n function createNdarray(style, buffer, matrix) {\n var shape = matrix.shape;\n var stride = matrix.stride;\n var offset = 0;\n return style === 'scijs'\n ? ndarray(buffer, shape, stride, offset)\n : stdlib.array(buffer, {\n shape,\n stride,\n offset,\n order: 'column-major'\n });\n }\n return function showMatrix(matrix, opts) {\n if (typeof opts === 'string') {\n opts = Object.assign(arguments[1] || {}, { label: opts });\n }\n var nd = createNdarray('scijs', matrix.buffer, matrix);\n var view = showNdarray(nd, opts);\n view.value = matrix;\n return view;\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2777,"value":"showMatrix(A)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2783,"value":"viewof blockAndDiag = md`## Block and diagonal operations`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3052,"value":"viewof diagonals = md`### Diagonals`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2892,"value":"md`You may use \\`.diagonal()\\` to reference a diagonal. This function returns a *view* of the band which you may use to modify values.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2786,"value":"Eigen.ensure(invalidation, () => {\n var A = Matrixf.zeros([8, 8]);\n\n // Set diagonal bands:\n A.diagonal(3).setConstant(-5);\n A.diagonal(0).setConstant(20);\n A.diagonal(-3).assign(Matrixf.columnVector([12, 13, 14, 15, 16]));\n\n return A.show();\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":3055,"value":"viewof blocks = md`### Blocks`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2895,"value":"md`You may use \\`block(startRow, startCol, shape)\\` to reference a submatrix. In a similar manner to \\`diagonal\\`, this does not copy the block but instead allows you to modify or mutate values directly.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2795,"value":"Eigen.ensure(invalidation, () => {\n var B = Matrixf.zeros([5, 5]);\n\n B.block(0, 1, [2, 2]).setConstant(25);\n\n B.block(1, 3, [2, 2]).assign(Matrixf.fromArray([10, 11, 12, 13], [2, 2]));\n\n B.bottomLeftCorner([1, 2]).setConstant(35);\n B.bottomRightCorner([2, 1]).setConstant(42);\n B.topRightCorner([1, 1]).setConstant(57);\n\n return B.show();\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":2923,"value":"md`Blocks are merely references and only have limited functionality. To use a block more generally, you may copy it into a new matrix using \\`Matrixf.fromBlock\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2865,"value":"Eigen.ensure(invalidation, () => {\n // Construct a large matrix:\n var B = Matrixf.random([20, 20]);\n\n // Define a block which references the original:\n var submatrix = B.block(5, 2, [3, 2]);\n\n // Construct and return a new matrix:\n return Matrixf.fromBlock(submatrix).show();\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":3059,"value":"viewof rowsCols = md`### Rows and columns`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3009,"value":"md`Rows and columns work similarly. You can \\`setZero()\\`, \\`setConstant()\\`, \\`setRandom()\\`, \\`set()\\` individual values, or \\`assign()\\` a vector.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3004,"value":"Eigen.ensure(invalidation, () => {\n // Construct a large matrix:\n var A = Matrixf.zeros([4, 4]);\n A.col(0).setConstant(5);\n A.row(2).setConstant(4);\n A.col(3).assign(Matrixf.columnVector([1, 2, 3, 4]));\n\n return A.show();\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":2885,"value":"md`You may use \\`.map()\\` to return a copy of a matrix with the function applied.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2874,"value":"A.map((row, col, value) => value + 2, invalidation).show()","pinned":true,"hidden":false,"mode":"ojs"},{"id":2909,"value":"viewof norms = md`## Norms`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2901,"value":"A.squaredNorm()","pinned":true,"hidden":false,"mode":"ojs"},{"id":2903,"value":"A.norm() // (square root of .squaredNorm())","pinned":true,"hidden":false,"mode":"ojs"},{"id":2905,"value":"A.l1Norm()","pinned":true,"hidden":false,"mode":"ojs"},{"id":2907,"value":"A.lInfNorm()","pinned":true,"hidden":false,"mode":"ojs"},{"id":2935,"value":"viewof advancedMemoryManagement = md`## Advanced memory management\n\nThe \\`Eigen.ensure(invalidation, () => ...)\\` block streamlines memory management, but there is one additional option worth mentioning. This construct frees memory allocated within _when the cell is invalidated_, but it may often be the case that the block allocates a number of temporary items but only outputs one.\n\nFor this case, \\`Eigen.ensure\\` passes as the first argument to the block an additional promise that resolves as soon as the block finishes executing. This may be used to immediately reclaim memory rather than leaving it around until the cell is invalidated at some time in the future.\n\nThe example below calls this immediately resolving promise \\`tmp\\` and uses it to immediately free all but the final result.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2943,"value":"Eigen.ensure(invalidation, tmp =>\n svd\n .matrixU(tmp)\n .mul(Matrixf.fromDiagonal(svd.singularValues(tmp), tmp), tmp)\n .mul(svd.matrixV(tmp).transpose(tmp), tmp)\n .map((i, j, value) => Math.round(value))\n).show('A = U \\\\Sigma V^T')","pinned":true,"hidden":false,"mode":"ojs"},{"id":2949,"value":"md`The behavior of this is subject to change. I'm debating whether immediate freeing is a better default and the special case should in fact be keeping the final result. This reduces typing but also makes things more error-prone since it would require action in order not to free the result as well. It's alternatively possible to check the \\`length\\` of the block to see if it has received the promise, but that starts to seem rather unexpected and convoluted.\n\nOf course it goes without saying that I wish you didn't have to manage memory at all.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":661,"value":"md`## Additional imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":271,"value":"import { slider, radio } from '@jashkenas/inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":636,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":2306,"value":"md`## License`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2310,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}