Menu
Search

Library Reference

overview

Report a bug
If you spot a problem with this page, click here to create a GitHub issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

kaleidic.lubeck

Lubeck - Linear Algebra

Authors:
Ilya Yaroshenko, Lars Tandle Kyllingstad (SciD author)
template BlasType(Iterators...)
Gets the type that can be used with Blas routines that all types can be implicitly converted to.
Slice!(BlasType!(IteratorA, IteratorB)*, 2) mtimes(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, 2, kindA) a, Slice!(IteratorB, 2, kindB) b);
General matrix-matrix multiplication. Allocates result to an uninitialized slice using GC.
Parameters:
Slice!(IteratorA, 2, kindA) a m(rows) x k(cols) matrix
Slice!(IteratorB, 2, kindB) b k(rows) x n(cols) matrix

Result m(rows) x n(cols)

Examples: 
import mir.ndslice;

auto a =
    [-5,  1,  7, 7, -4,
     -1, -5,  6, 3, -3,
     -5, -2, -3, 6,  0].sliced(3, 5);

auto b = slice!double(5, 4);
b[] =
    [[-5, -3,  3,  1],
     [ 4,  3,  6,  4],
     [-4, -2, -2,  2],
     [-1,  9,  4,  8],
     [ 9, 8,  3, -2]];

assert(mtimes(a, b) ==
    [[-42,  35,  -7, 77],
     [-69, -21, -42, 21],
     [ 23,  69,   3, 29]]
    );
Examples:
ger specialized case in mtimes 
import mir.ndslice;

// from https://github.com/kaleidicassociates/lubeck/issues/8
{
    auto a = [1.0f, 2.0f].sliced(2, 1);
    auto b = [1.0f, 2.0f].sliced(2, 1);
    assert(mtimes(a, b.transposed) == [[1, 2], [2, 4]]);
}
{
    auto a = [1.0, 2.0].sliced(1, 2);
    auto b = [1.0, 2.0].sliced(1, 2);
    assert(mtimes(a.transposed, b) == [[1, 2], [2, 4]]);
}
Examples: 
import mir.ndslice;

// from https://github.com/kaleidicassociates/lubeck/issues/3
Slice!(float*, 2) a = slice!float(1, 1);
Slice!(float*, 2, Universal) b1 = slice!float(16, 1).transposed;
Slice!(float*, 2) b2 = slice!float(1, 16);

a[] = 3;
b1[] = 4;
b2[] = 4;

// Confirm that this message does not appear
// Outputs: ** On entry to SGEMM  parameter number  8 had an illegal value
assert(a.mtimes(b1) == a.mtimes(b2));
Slice!(BlasType!(IteratorA, IteratorB)*) mtimes(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, 2, kindA) a, Slice!(IteratorB, 1, kindB) b);
General matrix-matrix multiplication. Allocates result to an uninitialized slice using GC.
Parameters:
Slice!(IteratorA, 2, kindA) a m(rows) x k(cols) matrix
Slice!(IteratorB, 1, kindB) b k(rows) x 1(cols) vector

Result m(rows) x 1(cols)

Slice!(BlasType!(IteratorA, IteratorB)*) mtimes(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorB, 1, kindB) a, Slice!(IteratorA, 2, kindA) b);
General matrix-matrix multiplication.
Parameters:
Slice!(IteratorB, 1, kindB) a 1(rows) x k(cols) vector
Slice!(IteratorA, 2, kindA) b k(rows) x n(cols) matrix

Result 1(rows) x n(cols)

Examples: 
import mir.ndslice;

auto a =
    [-5,  1,  7, 7, -4,
     -1, -5,  6, 3, -3,
     -5, -2, -3, 6,  0]
        .sliced(3, 5)
        .universal
        .transposed;

auto b = slice!double(5);
b[] = [-5, 4,-4,-1, 9];

assert(mtimes(b, a) == [-42, -69, 23]);
CommonType!(BlasType!IteratorA, BlasType!IteratorB) mtimes(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorB, 1, kindB) a, Slice!(IteratorA, 1, kindA) b);
Vector-vector multiplication (dot product).
Parameters:
Slice!(IteratorB, 1, kindB) a 1(rows) x k(cols) vector
Slice!(IteratorA, 1, kindA) b k(rows) x 1(cols) matrix

Result scalar

Examples: 
import mir.ndslice;

auto a = [1, 2, 4].sliced;
auto b = [3, 4, 2].sliced;
assert(a.mtimes(b) == 19);
auto inv(Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) a);
Calculates the inverse of a matrix.
Examples: 
import mir.complex;
import mir.ndslice;

auto a =  [
    1, 0, 2,
    2, 2, 0,
    0, 1, 1]
    .sliced(3, 3);

enum : double { _13 = 1.0/3.0, _16 = 1.0/6.0, _23 = 2.0/3.0 }
auto ans = [
    _13, _13, -_23,
    -_13,_16, _23,
    _13, -_16, _13]
    .sliced(a.shape);

import mir.algorithm.iteration: equal;
import mir.math.common: approxEqual;
assert(equal!((a, b) => a.approxEqual(b, 1e-10L, 1e-10L))(a.inv, ans));
assert(equal!((a, b) => a.approxEqual(b, 1e-10L, 1e-10L))(a.map!(a => Complex!double(a, 0)).inv.member!"re", ans));
Examples: 
import mir.ndslice.topology: iota;

try
{
    auto m = [3, 3].iota.inv;
    assert (false, "Matrix should be detected as singular");
}
catch (Exception e)
{
    assert (true);
}
struct SvdResult(T);
Slice!(T*, 2) u;
Slice!(T*) sigma;
Slice!(T*, 2) vt;
auto svd(Flag!"allowDestroy" allowDestroy = No.allowDestroy, string algorithm = "gesvd", SliceKind kind, Iterator)(Slice!(Iterator, 2, kind) matrix, Flag!"slim" slim = No.slim)
if (algorithm == "gesvd" || algorithm == "gesdd");
Computes the singular value decomposition.
Parameters:
Slice!(Iterator, 2, kind) matrix input M x N matrix
Flag!"slim" slim If true the first min(M,N) columns of u and the first min(M,N) rows of vt are returned in the ndslices u and vt.
Returns:
SvdResult. Results are allocated by the GC.
Examples: 
import mir.ndslice;

auto a =  [
     7.52,  -1.10,  -7.95,   1.08,
    -0.76,   0.62,   9.34,  -7.10,
     5.13,   6.62,  -5.66,   0.87,
    -4.75,   8.52,   5.75,   5.30,
     1.33,   4.91,  -5.49,  -3.52,
    -2.40,  -6.77,   2.34,   3.95]
    .sliced(6, 4);

auto r = a.svd;

auto sigma = slice!double(a.shape, 0);
sigma.diagonal[] = r.sigma;
auto m = r.u.mtimes(sigma).mtimes(r.vt);

import mir.algorithm.iteration: equal;
import mir.math.common: approxEqual;
assert(equal!((a, b) => a.approxEqual(b, 1e-8, 1e-8))(a, m));
Examples: 
import std.typecons: Yes;
import mir.ndslice;

auto a =  [
     7.52,  -1.10,  -7.95,   1.08,
    -0.76,   0.62,   9.34,  -7.10,
     5.13,   6.62,  -5.66,   0.87,
    -4.75,   8.52,   5.75,   5.30,
     1.33,   4.91,  -5.49,  -3.52,
    -2.40,  -6.77,   2.34,   3.95]
    .sliced(6, 4);

auto r = a.svd(Yes.slim);
assert(r.u.shape == [6, 4]);
assert(r.vt.shape == [4, 4]);
Slice!(BlasType!(IteratorA, IteratorB)*, 2) mldivide(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, 2, kindA) a, Slice!(IteratorB, 2, kindB) b);

Slice!(BlasType!(IteratorA, IteratorB)*) mldivide(IteratorA, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, 2, kindA) a, Slice!(IteratorB, 1, kindB) b);
Solve systems of linear equations AX = B for X. Computes minimum-norm solution to a linear least squares problem if A is not a square matrix.
Examples:
AX=B 
import mir.complex;
import std.meta: AliasSeq;
import mir.ndslice;

foreach(C; AliasSeq!(double, Complex!double))
{
    static if(is(C == Complex!double))
        alias transform = a => C(a, 0);
    else
        enum transform = "a";

    auto a = [
        1, -1,  1,
        2,  2, -4,
        -1,  5,  0].sliced(3, 3).map!transform;
    auto b = [
        2.0,  0,
        -6  , -6,
        9  ,  1].sliced(3, 2).map!transform;
    auto t = [
        1.0, -1,
        2  ,  0,
        3  ,  1].sliced(3, 2).map!transform;

    auto x = mldivide(a, b);
    assert(x == t);
}
Examples:
Ax=B 
import mir.complex;
import std.meta: AliasSeq;
import mir.ndslice;

foreach(C; AliasSeq!(double, Complex!double))
{
    static if(is(C == Complex!double))
        alias transform = a => C(a, 0);
    else
        enum transform = "a";

    auto a = [
        1, -1,  1,
        2,  2, -4,
        -1,  5,  0].sliced(3, 3).map!transform;
    auto b = [
        2.0,
        -6  ,
        9  ].sliced(3).map!transform;
    auto t = [
        1.0,
        2  ,
        3  ].sliced(3).map!transform;

    auto x = mldivide(a, b);
    assert(x == t);
}
Examples:
Least-Squares Solution of Underdetermined System 
import mir.complex;
import std.meta: AliasSeq;
import mir.ndslice;

foreach(C; AliasSeq!(double, )) //Complex!double fails for DMD>=2085
{
    static if(is(C == Complex!double))
        alias transform = a => C(a, 0);
    else
        enum transform = "a";

    auto a = [
        -0.57,  -1.28,  -0.39,   0.25,
        -1.93,   1.08,  -0.31,  -2.14,
        2.30,   0.24,   0.40,  -0.35,
        -1.93,   0.64,  -0.66,   0.08,
        0.15,   0.30,   0.15,  -2.13,
        -0.02,   1.03,  -1.43,   0.50,
        ].sliced(6, 4).map!transform;

    auto b = [
    -2.67,
    -0.55,
    3.34,
    -0.77,
    0.48,
    4.10,
    ].sliced.map!transform;

    auto x = [
        1.5339,
        1.8707,
        -1.5241,
        0.0392].sliced.map!transform;

    import mir.math.common: approxEqual;
    import mir.algorithm.iteration: all;
    alias appr = all!((a, b) => approxEqual(a, b, 1e-3, 1e-3));
    assert(appr(a.mldivide(b), x));
}
struct PcaResult(T);
Principal component analises result.
Slice!(T*, 2) coeff;
Principal component coefficients, also known as loadings.
Slice!(T*, 2) score;
Principal component scores.
Slice!(T*) latent;
Principal component variances.
auto pca(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) matrix, in Flag!"centerColumns" cc = Yes.centerColumns);
Principal component analysis of raw data.
Parameters:
Slice!(Iterator, 2, kind) matrix input M x N matrix, where 'M (rows)>= N(cols)'
Flag!"centerColumns" cc Flag to centern columns. True by default.
Returns:
PcaResult
Examples: 
import mir.ndslice;

import mir.math.common: approxEqual;
import mir.algorithm.iteration: equal;

auto ingedients = [
     7,  26,   6,  60,
     1,  29,  15,  52,
    11,  56,   8,  20,
    11,  31,   8,  47,
     7,  52,   6,  33,
    11,  55,   9,  22,
     3,  71,  17,   6,
     1,  31,  22,  44,
     2,  54,  18,  22,
    21,  47,   4,  26,
     1,  40,  23,  34,
    11,  66,   9,  12,
    10,  68,   8,  12].sliced(13, 4);

auto res = ingedients.pca;

auto coeff = [
    -0.067799985695474,  -0.646018286568728,   0.567314540990512,   0.506179559977705,
    -0.678516235418647,  -0.019993340484099,  -0.543969276583817,   0.493268092159297,
     0.029020832106229,   0.755309622491133,   0.403553469172668,   0.515567418476836,
     0.730873909451461,  -0.108480477171676,  -0.468397518388289,   0.484416225289198,
].sliced(4, 4);

auto score = [
     36.821825999449700,  -6.870878154227367,  -4.590944457629745,   0.396652582713912,
     29.607273420710964,   4.610881963526308,  -2.247578163663940,  -0.395843536696492,
    -12.981775719737618,  -4.204913183175938,   0.902243082694698,  -1.126100587210615,
     23.714725720918022,  -6.634052554708721,   1.854742000806314,  -0.378564808384691,
     -0.553191676624597,  -4.461732123178686,  -6.087412652325177,   0.142384896047281,
    -10.812490833309816,  -3.646571174544059,   0.912970791674604,  -0.134968810314680,
    -32.588166608817929,   8.979846284936063,  -1.606265913996588,   0.081763927599947,
     22.606395499005586,  10.725906457369449,   3.236537714483416,   0.324334774646368,
     -9.262587237675838,   8.985373347478788,  -0.016909578102172,  -0.543746175981799,
     -3.283969329640680, -14.157277337500918,   7.046512994833761,   0.340509860960606,
      9.220031117829379,  12.386080787220454,   3.428342878284624,   0.435152769664895,
    -25.584908517429557,  -2.781693148152386,  -0.386716066864491,   0.446817950545605,
    -26.903161834677597,  -2.930971165042989,  -2.445522630195304,   0.411607156409658,
].sliced(13, 4);

auto latent = [5.177968780739053, 0.674964360487231, 0.124054300480810, 0.002371532651878].sliced;
latent[] *= 100;

assert(equal!approxEqual(res.coeff, coeff));
assert(equal!approxEqual(res.score, score));
assert(equal!approxEqual(res.latent, latent));
Slice!(BlasType!Iterator*, 2) pinv(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) matrix, double tolerance = (double).nan);
Computes Moore-Penrose pseudoinverse of matrix.
Parameters:
Slice!(Iterator, 2, kind) matrix Input M x N matrix.
double tolerance The computation is based on SVD and any singular values less than tolerance are treated as zero.
Returns:
Moore-Penrose pseudoinverse matrix
Examples: 
import mir.ndslice;

auto a = [
    64,  2,  3, 61, 60,  6,
     9, 55, 54, 12, 13, 51,
    17, 47, 46, 20, 21, 43,
    40, 26, 27, 37, 36, 30,
    32, 34, 35, 29, 28, 38,
    41, 23, 22, 44, 45, 19,
    49, 15, 14, 52, 53, 11,
     8, 58, 59,  5,  4, 62].sliced(8, 6);

auto b = a.pinv;

auto result = [
    0.0177, -0.0165, -0.0164,  0.0174,  0.0173, -0.0161, -0.0160,  0.0170,
   -0.0121,  0.0132,  0.0130, -0.0114, -0.0112,  0.0124,  0.0122, -0.0106,
   -0.0055,  0.0064,  0.0060, -0.0043, -0.0040,  0.0049,  0.0045, -0.0028,
   -0.0020,  0.0039,  0.0046, -0.0038, -0.0044,  0.0064,  0.0070, -0.0063,
   -0.0086,  0.0108,  0.0115, -0.0109, -0.0117,  0.0139,  0.0147, -0.0141,
    0.0142, -0.0140, -0.0149,  0.0169,  0.0178, -0.0176, -0.0185,  0.0205].sliced(6, 8);

import mir.math.common: approxEqual;
import mir.algorithm.iteration: all;

assert(b.all!((a, b) => approxEqual(a, b, 1e-2, 1e-2))(result));
Slice!(BlasType!Iterator*, 2) cov(Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) matrix);
Covariance matrix.
Parameters:
Slice!(Iterator, 2, kind) matrix matrix whose rows represent observations and whose columns represent random variables.
Returns:
Normalized by N-1 covariance matrix.
Examples: 
import mir.ndslice;

import std.stdio;
import mir.ndslice;

auto c = 8.magic[0..$-1].cov;

auto result = [
     350.0000, -340.6667, -331.3333,  322.0000,  312.6667, -303.3333, -294.0000,  284.6667,
    -340.6667,  332.4762,  324.2857, -316.0952, -307.9048,  299.7143,  291.5238, -283.3333,
    -331.3333,  324.2857,  317.2381, -310.1905, -303.1429,  296.0952,  289.0476, -282.0000,
     322.0000, -316.0952, -310.1905,  304.2857,  298.3810, -292.4762, -286.5714,  280.6667,
     312.6667, -307.9048, -303.1429,  298.3810,  293.6190, -288.8571, -284.0952,  279.3333,
    -303.3333,  299.7143,  296.0952, -292.4762, -288.8571,  285.2381,  281.6190, -278.0000,
    -294.0000,  291.5238,  289.0476, -286.5714, -284.0952,  281.6190,  279.1429, -276.6667,
     284.6667, -283.3333, -282.0000,  280.6667,  279.3333, -278.0000, -276.6667,  275.3333].sliced(8, 8);
import mir.math.common: approxEqual;
import mir.algorithm.iteration: all;
assert(c.all!((a, b) => approxEqual(a, b, 1e-5, 1e-5))(result));
auto detSymmetric(Iterator, SliceKind kind)(char store, Slice!(Iterator, 2, kind) a);

auto det(Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) a);
Matrix determinant.
Examples: 
import mir.ndslice;
import mir.math;

// Check for zero-determinant shortcut.
auto ssing = [4, 2, 2, 1].sliced(2, 2);
auto ssingd = det(ssing);
assert (det(ssing) == 0);
assert (detSymmetric('L', ssing) == 0);

// check determinant of empty matrix
assert(slice!double(0, 0).det == 1);
// check determinant of zero matrix
assert(repeat(0, 9).sliced(3, 3).det == 0);

// General dense matrix.
int dn = 101;
auto d = uninitSlice!double(dn, dn);
foreach (k; 0 .. dn)
foreach (l; 0 .. dn)
    d[k,l] = 0.5 * (k == l ? (k + 1) * (k + 1) + 1 : 2 * (k + 1) * (l + 1));

auto dd = det(d);
import mir.math.common: approxEqual;
assert (approxEqual(dd, 3.539152633479803e289, double.epsilon.sqrt));

// Symmetric packed matrix
auto spa = [ 1.0, -2, 3, 4, 5, -6, -7, -8, -9, 10].sliced.stairs!"+"(4);
auto sp = [spa.length, spa.length].uninitSlice!double;
import mir.algorithm.iteration: each;
sp.stairs!"+".each!"a[] = b"(spa);
assert (detSymmetric('L', sp).approxEqual(5874.0, double.epsilon.sqrt));
assert (detSymmetric('U', sp.universal.transposed).approxEqual(5874.0, double.epsilon.sqrt));
struct EigSymmetricResult(T);
Eigenvalues and eigenvectors POD.
See Also:
eigSymmetric.
Slice!(T*) values;
Eigenvalues
Slice!(T*, 2) vectors;
Eigenvectors stored in rows
auto eigSymmetric(Flag!"computeVectors" cv = Yes.computeVectors, Iterator, SliceKind kind)(char store, Slice!(Iterator, 2, kind) a);
Eigenvalues and eigenvectors of symmetric matrix.
Returns:
EigSymmetricResult
Examples: 
import mir.ndslice;

import mir.ndslice.slice: sliced;
import mir.ndslice.topology: universal, map;
import mir.ndslice.dynamic: transposed;
import mir.math.common: approxEqual;
import mir.algorithm.iteration: all;

auto a = [
    1.0000, 0.5000, 0.3333, 0.2500,
    0.5000, 1.0000, 0.6667, 0.5000,
    0.3333, 0.6667, 1.0000, 0.7500,
    0.2500, 0.5000, 0.7500, 1.0000].sliced(4, 4);

auto eigr = eigSymmetric('L', a);

alias appr = all!((a, b) => approxEqual(a, b, 1e-3, 1e-3));

assert(appr(eigr.values, [0.2078,0.4078,0.8482,2.5362]));

auto test = [
     0.0693, -0.4422, -0.8105, 0.3778,
    -0.3618,  0.7420, -0.1877, 0.5322,
     0.7694,  0.0486,  0.3010, 0.5614,
    -0.5219, -0.5014,  0.4662, 0.5088].sliced(4, 4).transposed;


foreach (i; 0 .. 4)
    assert(appr(eigr.vectors[i], test[i]) || appr(eigr.vectors[i].map!"-a", test[i]));
struct LUResult(T);
LUResult consist lu factorization.
Slice!(T*, 2, Canonical) lut;
Matrix in witch lower triangular is L part of factorization (diagonal elements of L are not stored), upper triangular is U part of factorization.
Slice!(lapackint*) ipiv;
The pivot indices, for 1 <= i <= min(M,N), row i of the matrix was interchanged with row ipiv(i).
@property auto l();
L part of the factorization.
@property auto u();
U part of the factorization.
auto solve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, size_t N, SliceKind kind)(char trans, Slice!(Iterator, N, kind) b);
auto luDecomp(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) a);
Computes LU factorization of a general 'M x N' matrix 'A' using partial pivoting with row interchanges. The factorization has the form: \A = P * L * U Where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n), and U is upper triangular (upper trapezoidal if m < n).
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(Iterator, 2, kind) a input 'M x N' matrix for factorization.
Returns:
LUResalt
auto luSolve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, SliceKind kindB, size_t N, IteratorB, IteratorLU)(char trans, Slice!(IteratorLU, 2, Canonical) lut, Slice!(lapackint*) ipiv, Slice!(IteratorB, N, kindB) b);
Solves a system of linear equations \A * X = B, or \A**T * X = B with a general 'N x N' matrix 'A' using the LU factorization computed by luDecomp.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(IteratorLU, 2, Canonical) lut factorization of matrix 'A', A = P * L * U.
Slice!(lapackint*) ipiv the pivot indices from luDecomp.
Slice!(IteratorB, N, kindB) b the right hand side matrix B.
char trans specifies the form of the system of equations: 'N': A * X = B (No transpose) 'T': AT * X = B (Transpose) 'C': AT * X = B (Conjugate transpose = Transpose)
Returns:
Return solve of the system linear equations.
Examples: 
import mir.math;
import mir.ndslice;

auto A =
    [ 1,  4, -3,  5,  6,
     -2,  8,  5,  7,  8,
      3,  4,  7,  9,  1,
      2,  4,  6,  3,  2,
      6,  8,  3,  5,  2 ]
        .sliced(5, 5)
        .as!double.slice
        .canonical;

import mir.random.variable;
import mir.random.algorithm;
auto B = randomSlice(uniformVar(-100, 100), 5, 100);

auto LU = A.luDecomp();
auto X = LU.solve('N', B);

import mir.algorithm.iteration: equal;
assert(equal!((a, b) => fabs(a - b) < 1e-12)(mtimes(A, X), B));
struct LDLResult(T);
Consist LDL factorization;
char uplo;
uplo = 'U': Upper triangle is stored; 'L': lower triangle is stored.
Slice!(T*, 2, Canonical) matrix;
Matrix in witch lower triangular matrix is 'L' part of factorization, diagonal is 'D' part.
Slice!(lapackint*) ipiv;
The pivot indices. If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged and D(k, k) is a '1 x 1' diagonal block. If ipiv(k) = ipiv(k + 1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1, k:k+1) is a '2 x 2' diagonal block.
auto solve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, SliceKind kindB, size_t N, IteratorB)(Slice!(IteratorB, N, kindB) b);
Return solves a system of linear equations \A * X = B, using LDL factorization.
auto ldlDecomp(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(char uplo, Slice!(Iterator, 2, kind) a);
Computes the factorization of a real symmetric matrix A using the Bunch-Kaufman diagonal pivoting method. The for of the factorization is: \A = LDL**T Where L is product if permutation and unit lower triangular matrices, and D is symmetric and block diagonal with '1 x 1' and '2 x 2' diagonal blocks.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(Iterator, 2, kind) a input symmetric 'n x n' matrix for factorization.
char uplo 'U': Upper triangle is stored; 'L': lower triangle is stored. Returns:LDLResult
auto ldlSolve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, SliceKind kindB, size_t N, IteratorB, IteratorA)(char uplo, Slice!(IteratorA, 2, Canonical) a, Slice!(lapackint*) ipiv, Slice!(IteratorB, N, kindB) b);
Solves a system of linear equations \A * X = B with symmetric matrix 'A' using the factorization \A = U * D * UT, or \A = L * D * LT computed by ldlDecomp.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(IteratorA, 2, Canonical) a 'LD' or 'UD' matrix computed by ldlDecomp.
Slice!(lapackint*) ipiv details of the interchanges and the block structure of D as determined by ldlDecomp.
Slice!(IteratorB, N, kindB) b the right hand side matrix.
char uplo specifies whether the details of the factorization are stored as an upper or lower triangular matrix: 'U': Upper triangular, form is \A = U * D * UT; 'L': Lower triangular, form is \A = L * D * LT.
Returns:
The solution matrix.
Examples: 
import mir.ndslice;

auto A =
    [ 2.07,  3.87,  4.20, -1.15,
      3.87, -0.21,  1.87,  0.63,
      4.20,  1.87,  1.15,  2.06,
     -1.15,  0.63,  2.06, -1.81 ]
        .sliced(4, 4)
        .as!double.slice
        .canonical;

import mir.random.variable;
import mir.random.algorithm;
auto B = randomSlice(uniformVar(-100, 100), 4, 100);

auto LDL = ldlDecomp('L', A);
auto X = LDL.solve(B);

import mir.math.common: approxEqual;
import mir.algorithm.iteration: equal;
alias appr = equal!((a, b) => approxEqual(a, b, 1e-5, 1e-5));
assert(appr(mtimes(A, X), B));
auto choleskyDecomp(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(char uplo, Slice!(Iterator, 2, kind) a);
Computs Cholesky decomposition of symmetric positive definite matrix 'A'. The factorization has the form: \A = UT * U, if UPLO = 'U', or \A = L * LT, if UPLO = 'L' Where U is an upper triangular matrix and L is lower triangular.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(Iterator, 2, kind) a symmetric 'N x N' matrix.
char uplo if uplo is Upper, then upper triangle of A is stored, else lower.
Returns:
choleskyResult
Examples: 
import mir.algorithm.iteration: equal, eachUploPair;
import mir.ndslice;
import mir.random.algorithm;
import mir.random.variable;
import mir.math.common: approxEqual;

auto A =
       [ 25, double.nan, double.nan,
         15, 18,  double.nan,
         -5,  0, 11 ]
         .sliced(3, 3);

auto B = randomSlice(uniformVar(-100, 100), 3, 100);

auto C = choleskyDecomp('L', A);
auto X = C.solve(B);

A.eachUploPair!"a = b";
assert(equal!approxEqual(mtimes(A, X), B));
auto choleskySolve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, SliceKind kindB, size_t N, IteratorB, IteratorC)(char uplo, Slice!(IteratorC, 2, Canonical) c, Slice!(IteratorB, N, kindB) b);
Solves a system of linear equations A * X = B with a symmetric matrix A using the Cholesky factorization: \A = UT * U or \A = L * LT computed by choleskyDecomp.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(IteratorC, 2, Canonical) c the triangular factor 'U' or 'L' from the Cholesky factorization \A = UT * U or \A = L * LT, as computed by choleskyDecomp.
Slice!(IteratorB, N, kindB) b the right hand side matrix.
char uplo 'U': Upper triangle of A is stored; 'L': Lower triangle of A is stored.
Returns:
The solution matrix X.
Examples: 
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: as;
import std.typecons: Flag, Yes;

auto A =
        [ 1.0,  1,  3,
          1  ,  5,  5,
          3  ,  5, 19 ].sliced(3, 3);

auto B = [ 10.0,  157,  80 ].sliced;
auto C_ = B.dup.sliced(3, 1);

auto C = choleskyDecomp('U', A);
auto X = choleskySolve!(Yes.allowDestroy)(C.uplo, C.matrix, B);

import mir.math.common: approxEqual;
import mir.algorithm.iteration: equal;
alias appr = equal!((a, b) => approxEqual(a, b, 1e-5, 1e-5));
assert(appr(mtimes(A, X), C_));
struct QRResult(T);
Examples: 
import mir.ndslice;

auto A =
        [ 1,  1,  0,
          1,  0,  1,
          0,  1,  1 ]
          .sliced(3, 3)
          .as!double.slice;

auto Q_test =
        [ -0.7071068,  0.4082483,  -0.5773503,
          -0.7071068, -0.4082483,   0.5773503,
                   0,  0.8164966,   0.5773503]
          .sliced(3, 3)
          .as!double.slice;

auto R_test =
        [ -1.414214,  -0.7071068,   -0.7071068,
                  0,   1.2247449,    0.4082483,
                  0,           0,    1.1547005]
          .sliced(3, 3)
          .as!double.slice;

auto val = qrDecomp(A);

//saving these values to doublecheck they don't change later
auto val_matrix = val.matrix.slice;
auto val_tau = val.tau.slice;

import mir.math.common: approxEqual;
import mir.ndslice : equal;

auto r = val.R;
assert(equal!approxEqual(val.R, R_test));

auto q = val.Q;
assert(equal!approxEqual(val.Q, Q_test));

//double-checking values do not change
assert(equal!approxEqual(val_matrix, val.matrix));
assert(equal!approxEqual(val_tau, val.tau));

auto a = val.reconstruct;
assert(equal!approxEqual(A, a));
Slice!(T*, 2, Canonical) matrix;
Matrix in witch the elements on and above the diagonal of the array contain the min(M, N) x N upper trapezoidal matrix 'R' (R is upper triangular if m >= n). The elements below the diagonal, with the array tau, represent the orthogonal matrix 'Q' as product of min(m, n).
Slice!(T*) tau;
The scalar factors of the elementary reflectors
auto solve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) b);
Solve the least squares problem: \min ||A * X - B|| Using the QR factorization: \A = Q * R computed by qrDecomp.
auto Q(Flag!"allowDestroy" allowDestroy = No.allowDestroy);
Extract the Q matrix
auto R();
Extract the R matrix
auto reconstruct();
Reconstruct the original matrix given a QR decomposition
auto reconstruct(T, U)(T q, U r)
if (isMatrix!T && isMatrix!U);
Reconstruct the original matrix given a QR decomposition
auto qrDecomp(Flag!"allowDestroy" allowDestroy = No.allowDestroy, Iterator, SliceKind kind)(Slice!(Iterator, 2, kind) a);
Computes a QR factorization of matrix 'a'.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(Iterator, 2, kind) a initial matrix
Returns:
QRResult
auto qrSolve(Flag!"allowDestroy" allowDestroy = No.allowDestroy, SliceKind kindB, size_t N, IteratorB, IteratorA, IteratorT)(Slice!(IteratorA, 2, Canonical) a, Slice!IteratorT tau, Slice!(IteratorB, N, kindB) b);
Solve the least squares problem: \min ||A * X - B|| Using the QR factorization: \A = Q * R computed by qrDecomp.
Parameters:
allowDestroy flag to delete the source matrix.
Slice!(IteratorA, 2, Canonical) a detalis of the QR factorization of the original matrix as returned by qrDecomp.
Slice!IteratorT tau details of the orhtogonal matrix Q.
Slice!(IteratorB, N, kindB) b right hand side matrix.
Returns:
solution matrix.
Examples: 
import mir.ndslice;

auto A =
        [ 3,  1, -1,  2,
         -5,  1,  3, -4,
          2,  0,  1, -1,
          1, -5,  3, -3 ]
          .sliced(4, 4)
          .as!double.slice;

import mir.random.variable;
import mir.random.algorithm;
auto B = randomSlice(uniformVar(-100, 100), 4, 100);

auto C = qrDecomp(A);
auto X = C.solve(B);

import mir.math.common: approxEqual;
import mir.algorithm.iteration: equal;
alias appr = equal!((a, b) => approxEqual(a, b, 1e-5, 1e-5));
assert(appr(mtimes(A, X), B));