1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
#! /bin/perl
#
# From: Reuben Settergren <rjs@jhu.edu>
#
# I wrote a perl script to randomly generate matrices with smallish integers,
# and test for a suitable determinant (typically +/-1), and ouput the matrix
# and its inverse.
#
# There are switches to make it -s[ymmetric] and/or -d[iagonally strong],
# which are necessary and good properties for covariance matrices. You may
# need to grab Math::MatrixReal and Math::Factor::XS from cpan.
use Math::MatrixReal;
use Math::Factor::XS qw(prime_factors);
use Getopt::Std;
$opt{n} = 1;
getopts('hsdn:', \%opt);
$usage .= "easyinverse.pl [-hsd] [-n N]\n";
$usage .= " -h this message\n";
$usage .= " -s symmetric\n";
$usage .= " -d diagonally dominant\n";
$usage .= " -n N number to generate (DEF 1)\n";
if ($opt{h}) {
print $usage;
exit;
}
sub randint{ # but not 0
my $min = shift or -5;
my $max = shift or +5;
my $ansa = 0;
while ($ansa == 0) {
my $r = rand;
$r *= ($max - $min + 1);
my $i = int($r);
$ansa = $min + $i;
}
return $ansa;
}
$size = 4;
sub printmateq {
my $lbl = shift;
my $mat = shift;
my $round = shift;
if ($round) {
for my $r (1..$size) {
for my $c (1..$size) {
my $elt = $mat->element($r,$c);
my $pelt = abs($elt);
my $prnd = int($pelt * $round + 0.5) / $round;
my $rnd = $prnd * ($elt < 0 ? -1 : 1);
$mat->assign($r, $c, $rnd);
}}
}
my @rows;
for my $r (1..$size) {
my @row = map {$mat->element($r, $_)} (1..$size);
my $rowstr = '{' . (join ',', @row) . '}';
push @rows, $rowstr;
}
my $matstr = '{' . (join ',', @rows) . '}';
$matstr =~ s!\.?000+\d+!!g;
print "$lbl = $matstr;\n";
}
for (1..$opt{n}) { # requested number of matrices
$M = new Math::MatrixReal($size,$size);
while (1) {
for $r (1..$size) {
if ($opt{d}) { $M->assign($r,$r, randint(5,10)) }
else { $M->assign($r,$r, randint(-3, 3)) }
for $c ($r+1..$size) {
$i = randint(-5,5);
$M->assign($r, $c, $i);
$M->assign($c, $r, ($opt{s} ? $i : randint(-5,5)));
}
}
$det = int($M->det + 0.5);
next if $det == 0;
@factors = prime_factors(abs($det));
$all_good_factors = 1; # cross your fingers!
for $f (@factors) {
if ($f != 2 || $f != 5) {
$all_good_factors = 0;
last;
}
}
next unless $all_good_factors;
# test for positive semidefinite ~ all positive eigenvalues
$evals = $M->sym_eigenvalues;
$all_positive_evals = 1; # keep hope alive!
for $r (1..$size) {
if ($evals->element($r,1) <= 0) {
$all_positive_evals = 0;
last;
}
}
last if $all_positive_evals;
}
printmateq(".mat", $M);
printmateq(".inv", $M->inverse, 1000);
}
|
