Add specification for returning the least-squares solution to a linear matrix equation (linalg: lstsq) (#119)

* Stub spec

* Document keyword

* Update spec

* Add missing parenthesis

* Remove duplicate words

* Rename keyword argument and add support for setting tolerance to a float

* Update copy

* Reorder sentences

* Rename keyword argument

* Fix name

* Update copy

* Add support for providing an ordinate vector

* Update dtype requirements

* Update dtype requirements

* Update type annotation

* Update copy

* Move API to submodule

Author: kgryte

spec/API_specification/linear_algebra_functions.md

@@ -147,10 +147,35 @@ Computes the multiplicative inverse of a square matrix (or a stack of square mat
 
     -   an array containing the multiplicative inverses. The returned array must have a floating-point data type determined by {ref}`type-promotion` and must have the same shape as `x`.
 
-(function-lstsq)=
-### lstsq()
+(function-linalg-lstsq)=
+### linalg.lstsq(x1, x2, /, *, rtol=None)
 
-TODO
+Returns the least-squares solution to a linear matrix equation `Ax = b`.
+
+#### Parameters
+
+-   **x1**: _<array>_
+
+    -   coefficient array `A` having shape `(..., M, N)` and whose innermost two dimensions form `MxN` matrices. Should have a floating-point data type.
+
+-   **x2**: _<array>_
+
+    -   ordinate (or "dependent variable") array `b`. If `x2` has shape `(..., M)`, `x2` is equivalent to an array having shape `(..., M, 1)`, and `shape(x2)` must be compatible with `shape(x1)[:-1]` (see {ref}`broadcasting`). If `x2` has shape `(..., M, K)`, each column `k` defines a set of ordinate values for which to compute a solution, and `shape(x2)[:-1]` must be compatible with `shape(x1)[:-1]` (see {ref}`broadcasting`). Should have a floating-point data type.
+
+-   **rtol**: _Optional\[ Union\[ float, <array> ] ]_
+
+    -   relative tolerance for small singular values. Singular values less than or equal to `rtol * largest_singular_value` are set to zero. If a `float`, the value is equivalent to a zero-dimensional array having a data type determined by {ref}`type-promotion` (as applied to `x1` and `x2`) and must be broadcast against each matrix. If an `array`, must have a floating-point data type and must be compatible with `shape(x1)[:-2]` (see {ref}`broadcasting`). If `None`, the default value is `max(M, N) * eps`, where `eps` must be the machine epsilon associated with the floating-point data type determined by {ref}`type-promotion` (as applied to `x1` and `x2`). Default: `None`.
+
+#### Returns
+
+-   **out**: _Tuple\[ <array>, <array>, <array>, <array> ]_
+
+    -   a namedtuple `(x, residuals, rank, s)` whose
+    
+        -   first element must have the field name `x` and must be an array containing the least-squares solution for each `MxN` matrix in `x1`. The array containing the solutions must have shape `(N, K)` and must have a floating-point data type determined by {ref}`type-promotion`.
+        -   second element must have the field name `residuals` and must be an array containing the sum of squares residuals (i.e., the squared Euclidean 2-norm for each column in `b - Ax`). The array containing the residuals must have shape `(K,)` and must have a floating-point data type determined by {ref}`type-promotion`.
+        -   third element must have the field name `rank` and must be an array containing the effective rank of each `MxN` matrix. The array containing the ranks must have shape `shape(x1)[:-2]` and must have an integer data type.
+        -   fourth element must have the field name `s` and must be an array containing the singular values for each `MxN` matrix in `x1`. The array containing the singular values must have shape `(..., min(M, N))` and must have a floating-point data type determined by {ref}`type-promotion`.
 
 (function-matmul)=
 ### matmul()