Code Coverage for /src/SciPhp/LinAlg/CholeskyTrait.php

1. /src/SciPhp
2. LinAlg
3. CholeskyTrait.php

1	<?php
2
3	declare(strict_types=1);
4
5	namespace SciPhp\LinAlg;
6
7	use SciPhp\Exception\Message;
8	use SciPhp\NdArray;
9	use SciPhp\NumPhp as np;
10	use Webmozart\Assert\Assert;
11
12	/**
13	* Choleski decomposition
14	*/
15	trait CholeskyTrait
16	{
17	/**
18	* Decomposition of a Hermitian, positive-definite matrix into the
19	* product of a lower triangular matrix
20	*
21	* @param  \SciPhp\NdArray|array $matrix
22	*
23	* @link http://sciphp.org/linalg.cholesky
24	* Documentation for cholesky()
25	*
26	* @since 0.3.0
27	* @api
28	*/
29	final public function cholesky($matrix): NdArray
30	{
31	np::transform($matrix, true);
32
33	$shape = $matrix->shape;
34
35	Assert::true(
36	$matrix->is_square(),
37	sprintf(
38	Message::MAT_NOT_SQUARE,
39	trim((string) np::ar($shape))
40	)
41	);
42
43	$size = $shape[0];
44	$w = $matrix->copy();
45
46	$l = np::zeros($shape);
47
48	for ($a = 0; $a < $size; $a++) {
49	Assert::greaterThan(
50	$w[$a][$a],
51	0,
52	'Not a positive-definite matrix'
53	);
54
55	$l["{$a}, {$a}"] = sqrt($w[$a][$a]);
56
57	for ($b = $a + 1; $b < $size; $b++) {
58	Assert::eq(
59	$matrix["{$b}, {$a}"],
60	$matrix["{$a}, {$b}"],
61	'Not a symetric matrix'
62	);
63
64	$l["{$b}, {$a}"] = $w["{$b}, {$a}"] / $l["{$a}, {$a}"];
65
66	for ($c = $a + 1; $c <= $b; $c++) {
67	$w["${b}, ${c}"] = $w["{$b}, {$c}"]
68	- $l["{$b}, {$a}"]
69	* $l["{$c}, {$a}"];
70	}
71	}
72	}
73
74	return $l;
75	}
76	}