line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
|
2
|
|
|
|
|
|
|
# |
3
|
|
|
|
|
|
|
# GENERATED WITH PDL::PP! Don't modify! |
4
|
|
|
|
|
|
|
# |
5
|
|
|
|
|
|
|
package PDL::SVDLIBC; |
6
|
|
|
|
|
|
|
|
7
|
|
|
|
|
|
|
@EXPORT_OK = qw( PDL::PP _svdccsencode svdlas2a PDL::PP svdlas2 svdlas2aw PDL::PP svdlas2w svdlas2ad PDL::PP svdlas2d PDL::PP svdindexND svdindexNDt PDL::PP svdindexccs PDL::PP svderror ); |
8
|
|
|
|
|
|
|
%EXPORT_TAGS = (Func=>[@EXPORT_OK]); |
9
|
|
|
|
|
|
|
|
10
|
2
|
|
|
2
|
|
561269
|
use PDL::Core; |
|
2
|
|
|
|
|
9
|
|
|
2
|
|
|
|
|
11
|
|
11
|
2
|
|
|
2
|
|
408
|
use PDL::Exporter; |
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
10
|
|
12
|
2
|
|
|
2
|
|
48
|
use DynaLoader; |
|
2
|
|
|
|
|
3
|
|
|
2
|
|
|
|
|
155
|
|
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
$PDL::SVDLIBC::VERSION = 0.19; |
17
|
|
|
|
|
|
|
@ISA = ( 'PDL::Exporter','DynaLoader' ); |
18
|
|
|
|
|
|
|
push @PDL::Core::PP, __PACKAGE__; |
19
|
|
|
|
|
|
|
bootstrap PDL::SVDLIBC $VERSION; |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
=pod |
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
=head1 NAME |
27
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
PDL::SVDLIBC - PDL interface to Doug Rohde's SVD C Library |
29
|
|
|
|
|
|
|
|
30
|
|
|
|
|
|
|
=head1 SYNOPSIS |
31
|
|
|
|
|
|
|
|
32
|
|
|
|
|
|
|
use PDL; |
33
|
|
|
|
|
|
|
use PDL::SVDLIBC; |
34
|
|
|
|
|
|
|
|
35
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
36
|
|
|
|
|
|
|
## Input matrix (dense) |
37
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
38
|
|
|
|
|
|
|
$n = 100; ##-- number of columns |
39
|
|
|
|
|
|
|
$m = 50; ##-- number of rows |
40
|
|
|
|
|
|
|
$a = random(double,$n,$m); ##-- random matrix |
41
|
|
|
|
|
|
|
|
42
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
43
|
|
|
|
|
|
|
## Output pdls |
44
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
45
|
|
|
|
|
|
|
$d = $n; ##-- max number of output dimensions |
46
|
|
|
|
|
|
|
$ut = zeroes(double,$m,$d); ##-- left singular components |
47
|
|
|
|
|
|
|
$s = zeroes(double,$d); ##-- singular values (diagnonal vector) |
48
|
|
|
|
|
|
|
$vt = zeroes(double,$n,$d); ##-- right singular components |
49
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
51
|
|
|
|
|
|
|
## Singular Value Decomposition (dense) |
52
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
53
|
|
|
|
|
|
|
svdlas2d($a, $maxiters, $end, $kappa, $ut, $s, $vt); |
54
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
56
|
|
|
|
|
|
|
## Singular Value Decomposition (sparse, using direct whichND()-encoding) |
57
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
58
|
|
|
|
|
|
|
$which = whichND($a)->qsortvec(); |
59
|
|
|
|
|
|
|
$nzvals = indexND($a,$which); |
60
|
|
|
|
|
|
|
|
61
|
|
|
|
|
|
|
svdlas2w($which, $nzvals, $n, $m, $maxiters, $end, $kappa, $ut, $s, $vt); |
62
|
|
|
|
|
|
|
|
63
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
64
|
|
|
|
|
|
|
## Singular Value Decomposition (sparse, using PDL::CCS encoding) |
65
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
66
|
|
|
|
|
|
|
use PDL::CCS; |
67
|
|
|
|
|
|
|
($ptr,$rowids,$nzvals) = ccsencode($a); |
68
|
|
|
|
|
|
|
$ptr->reshape($ptr->nelem+1); |
69
|
|
|
|
|
|
|
$ptr->set(-1, $rowids->nelem); |
70
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
svdlas2($ptr, $rowids, $nzvals, $m, $maxiters, $end, $kappa, $ut, $s, $vt); |
72
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
74
|
|
|
|
|
|
|
## SVD decoding (lookup) |
75
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
76
|
|
|
|
|
|
|
$vals = svdindexND ($u, $s, $v, $which); |
77
|
|
|
|
|
|
|
$vals = svdindexNDt($ut,$s,$vt, $which); |
78
|
|
|
|
|
|
|
$vals = svdindexccs($u, $s, $v, $ptr,$rowids); |
79
|
|
|
|
|
|
|
$err = svderror ($u, $s, $v, $ptr,$rowids,$nzvals); |
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
=head1 DESCRIPTION |
82
|
|
|
|
|
|
|
|
83
|
|
|
|
|
|
|
PDL::SVDLIBC provides a PDL interface to the SVDLIBC routines |
84
|
|
|
|
|
|
|
for singular value decomposition of large sparse matrices. |
85
|
|
|
|
|
|
|
SVDLIBC is available from http://tedlab.mit.edu/~dr/SVDLIBC/ |
86
|
|
|
|
|
|
|
|
87
|
|
|
|
|
|
|
=cut |
88
|
|
|
|
|
|
|
|
89
|
|
|
|
|
|
|
|
90
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
|
93
|
|
|
|
|
|
|
|
94
|
|
|
|
|
|
|
|
95
|
|
|
|
|
|
|
=head1 FUNCTIONS |
96
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
|
98
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
=cut |
100
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
|
102
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
|
104
|
|
|
|
|
|
|
|
105
|
2
|
|
|
2
|
|
10
|
use strict; |
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
1206
|
|
106
|
|
|
|
|
|
|
|
107
|
|
|
|
|
|
|
=pod |
108
|
|
|
|
|
|
|
|
109
|
|
|
|
|
|
|
=head1 SVDLIBC Globals |
110
|
|
|
|
|
|
|
|
111
|
|
|
|
|
|
|
There are several global data structures still lurking in the |
112
|
|
|
|
|
|
|
SVDLIBC code, so expect problems if you are trying to run more |
113
|
|
|
|
|
|
|
than one 'las2' procedure at once (even in different processes). |
114
|
|
|
|
|
|
|
|
115
|
|
|
|
|
|
|
PDL::SVDLIBC provides access to (some of) the SVDLIBC globals |
116
|
|
|
|
|
|
|
through the following functions, which are not exported. |
117
|
|
|
|
|
|
|
|
118
|
|
|
|
|
|
|
=cut |
119
|
|
|
|
|
|
|
|
120
|
|
|
|
|
|
|
|
121
|
|
|
|
|
|
|
|
122
|
|
|
|
|
|
|
=pod |
123
|
|
|
|
|
|
|
|
124
|
|
|
|
|
|
|
=head2 PDL::SVDLIBC::verbosity() |
125
|
|
|
|
|
|
|
|
126
|
|
|
|
|
|
|
=head2 PDL::SVDLIBC::verbosity($level) |
127
|
|
|
|
|
|
|
|
128
|
|
|
|
|
|
|
Get/set the current SVDLIBC verbosity level. |
129
|
|
|
|
|
|
|
Valid values for $level are between 0 (no messages) and |
130
|
|
|
|
|
|
|
2 (many messages). |
131
|
|
|
|
|
|
|
|
132
|
|
|
|
|
|
|
=cut |
133
|
|
|
|
|
|
|
|
134
|
|
|
|
|
|
|
|
135
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
|
137
|
|
|
|
|
|
|
=pod |
138
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
=head2 PDL::SVDLIBC::svdVersion() |
140
|
|
|
|
|
|
|
|
141
|
|
|
|
|
|
|
Returns a string representing the SVDLIBC version |
142
|
|
|
|
|
|
|
this module was compiled with. |
143
|
|
|
|
|
|
|
|
144
|
|
|
|
|
|
|
=cut |
145
|
|
|
|
|
|
|
|
146
|
|
|
|
|
|
|
|
147
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
|
149
|
|
|
|
|
|
|
=pod |
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
=head1 SVD Utilities |
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
=cut |
154
|
|
|
|
|
|
|
|
155
|
|
|
|
|
|
|
|
156
|
|
|
|
|
|
|
|
157
|
|
|
|
|
|
|
|
158
|
|
|
|
|
|
|
|
159
|
|
|
|
|
|
|
=head2 _svdccsencode |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
=for sig |
162
|
|
|
|
|
|
|
|
163
|
|
|
|
|
|
|
Signature: (double a(n,m); indx [o]ptr(n1); indx [o]rowids(nnz); double [o]nzvals(nnz)) |
164
|
|
|
|
|
|
|
|
165
|
|
|
|
|
|
|
|
166
|
|
|
|
|
|
|
=for ref |
167
|
|
|
|
|
|
|
|
168
|
|
|
|
|
|
|
info not available |
169
|
|
|
|
|
|
|
|
170
|
|
|
|
|
|
|
|
171
|
|
|
|
|
|
|
=for bad |
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
_svdccsencode does not process bad values. |
174
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
175
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
|
177
|
|
|
|
|
|
|
=cut |
178
|
|
|
|
|
|
|
|
179
|
|
|
|
|
|
|
|
180
|
|
|
|
|
|
|
|
181
|
|
|
|
|
|
|
|
182
|
|
|
|
|
|
|
|
183
|
|
|
|
|
|
|
|
184
|
|
|
|
|
|
|
*_svdccsencode = \&PDL::_svdccsencode; |
185
|
|
|
|
|
|
|
|
186
|
|
|
|
|
|
|
|
187
|
|
|
|
|
|
|
|
188
|
|
|
|
|
|
|
|
189
|
|
|
|
|
|
|
=pod |
190
|
|
|
|
|
|
|
|
191
|
|
|
|
|
|
|
=head2 svdlas2a |
192
|
|
|
|
|
|
|
|
193
|
|
|
|
|
|
|
=for sig |
194
|
|
|
|
|
|
|
|
195
|
|
|
|
|
|
|
indx ptr(nplus1); |
196
|
|
|
|
|
|
|
indx rowids(nnz); |
197
|
|
|
|
|
|
|
double nzvals(nnz); |
198
|
|
|
|
|
|
|
indx m(); ##-- default: max($rowids)+1 |
199
|
|
|
|
|
|
|
int d(); ##-- default: max(nplus1-1,m) |
200
|
|
|
|
|
|
|
int iterations(); ##-- default: 2*$d |
201
|
|
|
|
|
|
|
double end(2); ##-- default: [-1e-30,1e-30] |
202
|
|
|
|
|
|
|
double kappa(); ##-- default: 1e-6 |
203
|
|
|
|
|
|
|
double [o]ut(m,d); ##-- default: new |
204
|
|
|
|
|
|
|
double [o] s(d); ##-- default: new |
205
|
|
|
|
|
|
|
double [o]vt(n,d); ##-- default: new |
206
|
|
|
|
|
|
|
|
207
|
|
|
|
|
|
|
Uses a variant of the single-vector Lanczos method (Lanczos, 1950) |
208
|
|
|
|
|
|
|
to compute the singular value decomposition of a sparse matrix with |
209
|
|
|
|
|
|
|
$m() rows and data encoded |
210
|
|
|
|
|
|
|
in Harwell-Boeing sparse format in the input parameters $ptr(), $rowids(), |
211
|
|
|
|
|
|
|
and $nzvals(). See L<"PDL::CCS"> for a way to acquire these parameters |
212
|
|
|
|
|
|
|
from a dense input matrix, but note that for svdlas2(), the |
213
|
|
|
|
|
|
|
column pointer $ptr() is of size ($n+1) for a dense matrix $a with |
214
|
|
|
|
|
|
|
$n columns, where $ptr($n)==$nnz is the total number of nonzero |
215
|
|
|
|
|
|
|
values in $a. |
216
|
|
|
|
|
|
|
|
217
|
|
|
|
|
|
|
$iterations() is the maximum number of Lanczos iterations to perform. |
218
|
|
|
|
|
|
|
|
219
|
|
|
|
|
|
|
$end() specifies two endpoints of an interval within which all unwanted |
220
|
|
|
|
|
|
|
eigenvalues lie. |
221
|
|
|
|
|
|
|
|
222
|
|
|
|
|
|
|
$kappa() is a double containing the relative accuracy of Ritz |
223
|
|
|
|
|
|
|
values acceptable as eigenvalues. |
224
|
|
|
|
|
|
|
|
225
|
|
|
|
|
|
|
The left singular components are returned in the matrix $ut(), |
226
|
|
|
|
|
|
|
the singular values themselved in the vector $s(), and the |
227
|
|
|
|
|
|
|
right singular components in the matrix $vt(). Note that |
228
|
|
|
|
|
|
|
$ut() and $vt() are transposed, and must be specified explicitly |
229
|
|
|
|
|
|
|
in the call, so that the degree of reduction (the size parameter $d) |
230
|
|
|
|
|
|
|
can be determined. If $d==$n, then a full decomposition |
231
|
|
|
|
|
|
|
will be computed, and on return, $ut() and $vt() should be transposed |
232
|
|
|
|
|
|
|
instances of the matrices $u() and $v() as returned by PDL::MatrixOps::svd(). |
233
|
|
|
|
|
|
|
|
234
|
|
|
|
|
|
|
The Lanczos method as used here seems to be consistently the |
235
|
|
|
|
|
|
|
fastest. This algorithm has the drawback that the low order singular |
236
|
|
|
|
|
|
|
values may be relatively imprecise, but that is not a problem for most |
237
|
|
|
|
|
|
|
users who only want the higher-order values or who can tolerate some |
238
|
|
|
|
|
|
|
imprecision. |
239
|
|
|
|
|
|
|
|
240
|
|
|
|
|
|
|
See also: svdlas2aw(), svdlas2d() |
241
|
|
|
|
|
|
|
|
242
|
|
|
|
|
|
|
=cut |
243
|
|
|
|
|
|
|
|
244
|
|
|
|
|
|
|
## ($iters,$end,$kappa,$ut,$s,$vt) = svddefaults($n=$nrows,$m=$ncols,$d, $iters,...) |
245
|
|
|
|
|
|
|
## + returns default values |
246
|
|
|
|
|
|
|
## + changed calling conventions in v0.14 |
247
|
|
|
|
|
|
|
## - WAS: svddefaults($nrows,$cols, $d,$iters,...) ##-- SVDLIBC-style (col-primary) |
248
|
|
|
|
|
|
|
## ~ svddefaults($m, $n, $d,$iters,...) ##-- SVDLIBC-style (for dense $a(n,m)) |
249
|
|
|
|
|
|
|
## - NOW: svddefaults($n, $m, $d,$iters,...) ##-- pdl-style |
250
|
|
|
|
|
|
|
sub svddefaults { |
251
|
7
|
|
|
7
|
0
|
19
|
my ($n,$m,$d, $iters,$end,$kappa,$ut,$s,$vt) = @_; |
252
|
7
|
50
|
|
|
|
23
|
$n = $n->at(0) if (UNIVERSAL::isa($n,'PDL')); |
253
|
7
|
50
|
|
|
|
23
|
$m = $m->at(0) if (UNIVERSAL::isa($m,'PDL')); |
254
|
7
|
50
|
|
|
|
19
|
$d = ($n >= $m ? $n : $m) if (!defined($d)); |
|
|
100
|
|
|
|
|
|
255
|
7
|
50
|
|
|
|
20
|
$iters = 2*$d if (!defined($iters)); |
256
|
7
|
50
|
|
|
|
25
|
$end = pdl(double,[-1e-30,1e-30]) if (!defined($end)); |
257
|
7
|
50
|
|
|
|
184
|
$kappa = pdl(double,1e-6) if (!defined($kappa)); |
258
|
7
|
50
|
|
|
|
371
|
$ut = PDL->zeroes(double,$m,$d) if (!defined($ut)); |
259
|
7
|
50
|
|
|
|
416
|
$s = PDL->zeroes(double, $d) if (!defined($s)); |
260
|
7
|
50
|
|
|
|
290
|
$vt = PDL->zeroes(double,$n,$d) if (!defined($vt)); |
261
|
7
|
|
|
|
|
311
|
return ($iters,$end,$kappa,$ut,$s,$vt); |
262
|
|
|
|
|
|
|
} |
263
|
|
|
|
|
|
|
|
264
|
|
|
|
|
|
|
sub svdlas2a { |
265
|
2
|
|
|
2
|
1
|
113795
|
my ($ptr,$rowids,$nzvals, $m,$d, @args) = @_; |
266
|
2
|
100
|
|
|
|
12
|
$m = $rowids->flat->max+1 if (!defined($m)); |
267
|
2
|
|
|
|
|
87
|
@args = svddefaults($ptr->dim(0)-1,$m,$d,@args); |
268
|
2
|
|
|
|
|
128
|
svdlas2($ptr,$rowids,$nzvals,$m,@args); |
269
|
2
|
|
|
|
|
22
|
return @args[3..5]; |
270
|
|
|
|
|
|
|
} |
271
|
|
|
|
|
|
|
|
272
|
|
|
|
|
|
|
|
273
|
|
|
|
|
|
|
|
274
|
|
|
|
|
|
|
|
275
|
|
|
|
|
|
|
|
276
|
|
|
|
|
|
|
=head2 svdlas2 |
277
|
|
|
|
|
|
|
|
278
|
|
|
|
|
|
|
=for sig |
279
|
|
|
|
|
|
|
|
280
|
|
|
|
|
|
|
Signature: ( |
281
|
|
|
|
|
|
|
indx ptr(nplus1); |
282
|
|
|
|
|
|
|
indx rowids(nnz); |
283
|
|
|
|
|
|
|
double nzvals(nnz); |
284
|
|
|
|
|
|
|
indx m(); |
285
|
|
|
|
|
|
|
int iterations(); |
286
|
|
|
|
|
|
|
double end(2); |
287
|
|
|
|
|
|
|
double kappa(); |
288
|
|
|
|
|
|
|
double [o]ut(m,d); |
289
|
|
|
|
|
|
|
double [o] s(d); |
290
|
|
|
|
|
|
|
double [o]vt(n,d); |
291
|
|
|
|
|
|
|
) |
292
|
|
|
|
|
|
|
|
293
|
|
|
|
|
|
|
|
294
|
|
|
|
|
|
|
Guts for svdlas2a(). |
295
|
|
|
|
|
|
|
No default instantiation, and slightly different calling conventions. |
296
|
|
|
|
|
|
|
|
297
|
|
|
|
|
|
|
|
298
|
|
|
|
|
|
|
=for bad |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
svdlas2 does not process bad values. |
301
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
302
|
|
|
|
|
|
|
|
303
|
|
|
|
|
|
|
|
304
|
|
|
|
|
|
|
=cut |
305
|
|
|
|
|
|
|
|
306
|
|
|
|
|
|
|
|
307
|
|
|
|
|
|
|
|
308
|
|
|
|
|
|
|
|
309
|
|
|
|
|
|
|
|
310
|
|
|
|
|
|
|
|
311
|
|
|
|
|
|
|
*svdlas2 = \&PDL::svdlas2; |
312
|
|
|
|
|
|
|
|
313
|
|
|
|
|
|
|
|
314
|
|
|
|
|
|
|
|
315
|
|
|
|
|
|
|
|
316
|
|
|
|
|
|
|
=pod |
317
|
|
|
|
|
|
|
|
318
|
|
|
|
|
|
|
=head2 svdlas2aw |
319
|
|
|
|
|
|
|
|
320
|
|
|
|
|
|
|
=for sig |
321
|
|
|
|
|
|
|
|
322
|
|
|
|
|
|
|
indx which(nnz,2); ##-- sorted indices of non-zero values |
323
|
|
|
|
|
|
|
double nzvals(nnz); ##-- non-zero values |
324
|
|
|
|
|
|
|
indx n(); ##-- default: max($indx(0,:))+1 |
325
|
|
|
|
|
|
|
indx m(); ##-- default: max($indx(1,:))+1 |
326
|
|
|
|
|
|
|
int d(); ##-- default: max(n,m) |
327
|
|
|
|
|
|
|
int iterations(); ##-- default: 2*$d |
328
|
|
|
|
|
|
|
double end(2); ##-- default: [-1e-30,1e-30] |
329
|
|
|
|
|
|
|
double kappa(); ##-- default: 1e-6 |
330
|
|
|
|
|
|
|
double [o]ut(m,d); ##-- default: new |
331
|
|
|
|
|
|
|
double [o] s(d); ##-- default: new |
332
|
|
|
|
|
|
|
double [o]vt(n,d); ##-- default: new |
333
|
|
|
|
|
|
|
|
334
|
|
|
|
|
|
|
As for svdlas2a(), but implicitly converts the index-encoded matrix |
335
|
|
|
|
|
|
|
($which(),$nzvals()) to an internal CCS-like sparse format |
336
|
|
|
|
|
|
|
before computing the decomposition. |
337
|
|
|
|
|
|
|
Should be slightly more efficient than using PDL::CCS::ccsencode() |
338
|
|
|
|
|
|
|
or similar if you already have $which() and $nzvals() available. |
339
|
|
|
|
|
|
|
These can be attained for a dense matric $a() e.g. by: |
340
|
|
|
|
|
|
|
|
341
|
|
|
|
|
|
|
$which = $a->whichND->qsortvec->xchg(0,1); |
342
|
|
|
|
|
|
|
$nzvals = $a->indexND($which->xchg(0,1)); |
343
|
|
|
|
|
|
|
|
344
|
|
|
|
|
|
|
For convenience, $which() will be implicitly transposed if it is passed |
345
|
|
|
|
|
|
|
as a list-of-vectors C<$whichND(2,nnz)> such as returned by L, |
346
|
|
|
|
|
|
|
but it must still be lexicographically sorted. |
347
|
|
|
|
|
|
|
|
348
|
|
|
|
|
|
|
See also: svdlas2a(), svdlas2d() |
349
|
|
|
|
|
|
|
|
350
|
|
|
|
|
|
|
=cut |
351
|
|
|
|
|
|
|
|
352
|
|
|
|
|
|
|
sub svdlas2aw { |
353
|
3
|
|
|
3
|
1
|
5898
|
my ($which,$nzvals, $n,$m,$d, @args) = @_; |
354
|
3
|
100
|
|
|
|
22
|
$which = $which->xchg(0,1) if ($which->dim(1) > $which->dim(0)); |
355
|
3
|
100
|
|
|
|
14
|
$n = $which->slice(":,0")->max+1 if (!defined($n)); |
356
|
3
|
100
|
|
|
|
135
|
$m = $which->slice(":,1")->max+1 if (!defined($m)); |
357
|
3
|
|
|
|
|
92
|
@args = svddefaults($n,$m,$d,@args); |
358
|
3
|
|
|
|
|
225
|
svdlas2w($which,$nzvals,$n,$m,@args); |
359
|
3
|
|
|
|
|
34
|
return @args[3..5]; |
360
|
|
|
|
|
|
|
} |
361
|
|
|
|
|
|
|
|
362
|
|
|
|
|
|
|
|
363
|
|
|
|
|
|
|
|
364
|
|
|
|
|
|
|
|
365
|
|
|
|
|
|
|
|
366
|
|
|
|
|
|
|
=head2 svdlas2w |
367
|
|
|
|
|
|
|
|
368
|
|
|
|
|
|
|
=for sig |
369
|
|
|
|
|
|
|
|
370
|
|
|
|
|
|
|
Signature: ( |
371
|
|
|
|
|
|
|
indx whichi(nnz,Two); |
372
|
|
|
|
|
|
|
double nzvals(nnz); |
373
|
|
|
|
|
|
|
indx n(); |
374
|
|
|
|
|
|
|
indx m(); |
375
|
|
|
|
|
|
|
int iterations(); |
376
|
|
|
|
|
|
|
double end(2); |
377
|
|
|
|
|
|
|
double kappa(); |
378
|
|
|
|
|
|
|
double [o]ut(m,d); |
379
|
|
|
|
|
|
|
double [o] s(d); |
380
|
|
|
|
|
|
|
double [o]vt(n,d); |
381
|
|
|
|
|
|
|
) |
382
|
|
|
|
|
|
|
|
383
|
|
|
|
|
|
|
|
384
|
|
|
|
|
|
|
Guts for svdlas2a(). |
385
|
|
|
|
|
|
|
No default instantiation, and slightly different calling conventions. |
386
|
|
|
|
|
|
|
|
387
|
|
|
|
|
|
|
|
388
|
|
|
|
|
|
|
=for bad |
389
|
|
|
|
|
|
|
|
390
|
|
|
|
|
|
|
svdlas2w does not process bad values. |
391
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
392
|
|
|
|
|
|
|
|
393
|
|
|
|
|
|
|
|
394
|
|
|
|
|
|
|
=cut |
395
|
|
|
|
|
|
|
|
396
|
|
|
|
|
|
|
|
397
|
|
|
|
|
|
|
|
398
|
|
|
|
|
|
|
|
399
|
|
|
|
|
|
|
|
400
|
|
|
|
|
|
|
|
401
|
|
|
|
|
|
|
*svdlas2w = \&PDL::svdlas2w; |
402
|
|
|
|
|
|
|
|
403
|
|
|
|
|
|
|
|
404
|
|
|
|
|
|
|
|
405
|
|
|
|
|
|
|
|
406
|
|
|
|
|
|
|
=pod |
407
|
|
|
|
|
|
|
|
408
|
|
|
|
|
|
|
=head2 svdlas2ad |
409
|
|
|
|
|
|
|
|
410
|
|
|
|
|
|
|
=for sig |
411
|
|
|
|
|
|
|
|
412
|
|
|
|
|
|
|
double a(n,m); |
413
|
|
|
|
|
|
|
int d(); ##-- default: max($n,$m) |
414
|
|
|
|
|
|
|
int iterations(); ##-- default: 2*$d |
415
|
|
|
|
|
|
|
double end(2); ##-- default: [-1e-30,1e-30] |
416
|
|
|
|
|
|
|
double kappa(); ##-- default: 1e-6 |
417
|
|
|
|
|
|
|
double [o]ut(m,d); ##-- default: new |
418
|
|
|
|
|
|
|
double [o] s(d); ##-- default: new |
419
|
|
|
|
|
|
|
double [o]vt(n,d); ##-- default: new |
420
|
|
|
|
|
|
|
|
421
|
|
|
|
|
|
|
As for svdlas2(), but implicitly converts the dense input matrix |
422
|
|
|
|
|
|
|
$a() to sparse format before computing the decomposition. |
423
|
|
|
|
|
|
|
|
424
|
|
|
|
|
|
|
=cut |
425
|
|
|
|
|
|
|
|
426
|
|
|
|
|
|
|
sub svdlas2ad { |
427
|
2
|
|
|
2
|
1
|
3534
|
my ($a,$d, @args) = @_; |
428
|
2
|
|
|
|
|
12
|
@args = svddefaults($a->dim(0),$a->dim(1),$d,@args); |
429
|
2
|
|
|
|
|
117
|
svdlas2d($a,@args); |
430
|
2
|
|
|
|
|
17
|
return @args[3..5]; |
431
|
|
|
|
|
|
|
} |
432
|
|
|
|
|
|
|
|
433
|
|
|
|
|
|
|
|
434
|
|
|
|
|
|
|
|
435
|
|
|
|
|
|
|
|
436
|
|
|
|
|
|
|
|
437
|
|
|
|
|
|
|
=head2 svdlas2d |
438
|
|
|
|
|
|
|
|
439
|
|
|
|
|
|
|
=for sig |
440
|
|
|
|
|
|
|
|
441
|
|
|
|
|
|
|
Signature: ( |
442
|
|
|
|
|
|
|
double a(n,m); |
443
|
|
|
|
|
|
|
int iterations(); |
444
|
|
|
|
|
|
|
double end(2); |
445
|
|
|
|
|
|
|
double kappa(); |
446
|
|
|
|
|
|
|
double [o]ut(m,d); |
447
|
|
|
|
|
|
|
double [o] s(d); |
448
|
|
|
|
|
|
|
double [o]vt(n,d); |
449
|
|
|
|
|
|
|
) |
450
|
|
|
|
|
|
|
|
451
|
|
|
|
|
|
|
|
452
|
|
|
|
|
|
|
Guts for _svdlas2d(). |
453
|
|
|
|
|
|
|
|
454
|
|
|
|
|
|
|
|
455
|
|
|
|
|
|
|
=for bad |
456
|
|
|
|
|
|
|
|
457
|
|
|
|
|
|
|
svdlas2d does not process bad values. |
458
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
459
|
|
|
|
|
|
|
|
460
|
|
|
|
|
|
|
|
461
|
|
|
|
|
|
|
=cut |
462
|
|
|
|
|
|
|
|
463
|
|
|
|
|
|
|
|
464
|
|
|
|
|
|
|
|
465
|
|
|
|
|
|
|
|
466
|
|
|
|
|
|
|
|
467
|
|
|
|
|
|
|
|
468
|
|
|
|
|
|
|
*svdlas2d = \&PDL::svdlas2d; |
469
|
|
|
|
|
|
|
|
470
|
|
|
|
|
|
|
|
471
|
|
|
|
|
|
|
|
472
|
|
|
|
|
|
|
|
473
|
|
|
|
|
|
|
|
474
|
|
|
|
|
|
|
=head2 svdindexND |
475
|
|
|
|
|
|
|
|
476
|
|
|
|
|
|
|
=for sig |
477
|
|
|
|
|
|
|
|
478
|
|
|
|
|
|
|
Signature: ( |
479
|
|
|
|
|
|
|
u(d,m); |
480
|
|
|
|
|
|
|
s(d); |
481
|
|
|
|
|
|
|
v(d,n); |
482
|
|
|
|
|
|
|
indx which(Two,nnz); |
483
|
|
|
|
|
|
|
[o] vals(nnz); |
484
|
|
|
|
|
|
|
) |
485
|
|
|
|
|
|
|
|
486
|
|
|
|
|
|
|
|
487
|
|
|
|
|
|
|
Lookup selected values in an SVD-encoded matrix, L-style. |
488
|
|
|
|
|
|
|
Should be equivalent to: |
489
|
|
|
|
|
|
|
|
490
|
|
|
|
|
|
|
($u x stretcher($s) x $v->xchg(0,1))->indexND($which) |
491
|
|
|
|
|
|
|
|
492
|
|
|
|
|
|
|
or its PDL-friendlier variant: |
493
|
|
|
|
|
|
|
|
494
|
|
|
|
|
|
|
($u * $s)->matmult($v->xchg(0,1))->indexND($which) |
495
|
|
|
|
|
|
|
|
496
|
|
|
|
|
|
|
... but only computes the specified values $which(), avoiding |
497
|
|
|
|
|
|
|
memory bottlenecks for large sparse matrices. |
498
|
|
|
|
|
|
|
This is a pure PDL::PP method, so you can use e.g. |
499
|
|
|
|
|
|
|
C for the SVD-encoded matrix if you wish. |
500
|
|
|
|
|
|
|
|
501
|
|
|
|
|
|
|
|
502
|
|
|
|
|
|
|
|
503
|
|
|
|
|
|
|
=for bad |
504
|
|
|
|
|
|
|
|
505
|
|
|
|
|
|
|
svdindexND does not process bad values. |
506
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
507
|
|
|
|
|
|
|
|
508
|
|
|
|
|
|
|
|
509
|
|
|
|
|
|
|
=cut |
510
|
|
|
|
|
|
|
|
511
|
|
|
|
|
|
|
|
512
|
|
|
|
|
|
|
|
513
|
|
|
|
|
|
|
|
514
|
|
|
|
|
|
|
|
515
|
|
|
|
|
|
|
|
516
|
|
|
|
|
|
|
*svdindexND = \&PDL::svdindexND; |
517
|
|
|
|
|
|
|
|
518
|
|
|
|
|
|
|
|
519
|
|
|
|
|
|
|
|
520
|
|
|
|
|
|
|
|
521
|
|
|
|
|
|
|
=pod |
522
|
|
|
|
|
|
|
|
523
|
|
|
|
|
|
|
=head2 svdindexNDt |
524
|
|
|
|
|
|
|
|
525
|
|
|
|
|
|
|
=for sig |
526
|
|
|
|
|
|
|
|
527
|
|
|
|
|
|
|
ut(m,d); s(d); vt(n,d); indx which(Two,nnz); [o] vals(nnz); |
528
|
|
|
|
|
|
|
|
529
|
|
|
|
|
|
|
Wrapper for L accepting transposed singular components |
530
|
|
|
|
|
|
|
$ut() and $vt() as returned by e.g. L. |
531
|
|
|
|
|
|
|
|
532
|
|
|
|
|
|
|
=cut |
533
|
|
|
|
|
|
|
|
534
|
|
|
|
|
|
|
sub svdindexNDt { |
535
|
2
|
|
|
2
|
1
|
2720
|
return svdindexND($_[0]->xchg(0,1),$_[1],$_[2]->xchg(0,1),@_[3..$#_]); |
536
|
|
|
|
|
|
|
} |
537
|
|
|
|
|
|
|
|
538
|
|
|
|
|
|
|
|
539
|
|
|
|
|
|
|
|
540
|
|
|
|
|
|
|
|
541
|
|
|
|
|
|
|
|
542
|
|
|
|
|
|
|
=head2 svdindexccs |
543
|
|
|
|
|
|
|
|
544
|
|
|
|
|
|
|
=for sig |
545
|
|
|
|
|
|
|
|
546
|
|
|
|
|
|
|
Signature: ( |
547
|
|
|
|
|
|
|
u(d,m); |
548
|
|
|
|
|
|
|
s(d); |
549
|
|
|
|
|
|
|
v(d,n); |
550
|
|
|
|
|
|
|
indx ptr(nplus1); |
551
|
|
|
|
|
|
|
indx rowids(nnz); |
552
|
|
|
|
|
|
|
[o] vals(nnz); |
553
|
|
|
|
|
|
|
) |
554
|
|
|
|
|
|
|
|
555
|
|
|
|
|
|
|
|
556
|
|
|
|
|
|
|
Lookup selected values in an SVD-encoded matrix using L-style indexing |
557
|
|
|
|
|
|
|
as for L. |
558
|
|
|
|
|
|
|
|
559
|
|
|
|
|
|
|
|
560
|
|
|
|
|
|
|
|
561
|
|
|
|
|
|
|
=for bad |
562
|
|
|
|
|
|
|
|
563
|
|
|
|
|
|
|
svdindexccs does not process bad values. |
564
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
565
|
|
|
|
|
|
|
|
566
|
|
|
|
|
|
|
|
567
|
|
|
|
|
|
|
=cut |
568
|
|
|
|
|
|
|
|
569
|
|
|
|
|
|
|
|
570
|
|
|
|
|
|
|
|
571
|
|
|
|
|
|
|
|
572
|
|
|
|
|
|
|
|
573
|
|
|
|
|
|
|
|
574
|
|
|
|
|
|
|
*svdindexccs = \&PDL::svdindexccs; |
575
|
|
|
|
|
|
|
|
576
|
|
|
|
|
|
|
|
577
|
|
|
|
|
|
|
|
578
|
|
|
|
|
|
|
|
579
|
|
|
|
|
|
|
|
580
|
|
|
|
|
|
|
=head2 svderror |
581
|
|
|
|
|
|
|
|
582
|
|
|
|
|
|
|
=for sig |
583
|
|
|
|
|
|
|
|
584
|
|
|
|
|
|
|
Signature: ( |
585
|
|
|
|
|
|
|
u(d,m); |
586
|
|
|
|
|
|
|
s(d); |
587
|
|
|
|
|
|
|
v(d,n); |
588
|
|
|
|
|
|
|
indx ptr(nplus1); |
589
|
|
|
|
|
|
|
indx rowids(nnz); |
590
|
|
|
|
|
|
|
nzvals(nnz); |
591
|
|
|
|
|
|
|
[o]err(); |
592
|
|
|
|
|
|
|
) |
593
|
|
|
|
|
|
|
|
594
|
|
|
|
|
|
|
|
595
|
|
|
|
|
|
|
Compute sum of squared errors for a sparse SVD-encoded matrix. |
596
|
|
|
|
|
|
|
Should be equivalent to: |
597
|
|
|
|
|
|
|
|
598
|
|
|
|
|
|
|
sum( ($a - ($u x stretcher($s) x $v->xchg(0,1)))**2 ) |
599
|
|
|
|
|
|
|
|
600
|
|
|
|
|
|
|
... but computes all values on-the-fly, avoiding |
601
|
|
|
|
|
|
|
memory bottlenecks for large sparse matrices. |
602
|
|
|
|
|
|
|
This is a pure PDL::PP method, so you can use e.g. |
603
|
|
|
|
|
|
|
C for the SVD-encoded matrix if you wish. |
604
|
|
|
|
|
|
|
|
605
|
|
|
|
|
|
|
Error contributions are computed even for "missing" (zero) values, |
606
|
|
|
|
|
|
|
so running time is O(n*m). |
607
|
|
|
|
|
|
|
Consider using L or L |
608
|
|
|
|
|
|
|
to compute error rates |
609
|
|
|
|
|
|
|
only for non-missing values if you have a large sparse matrix, e.g.: |
610
|
|
|
|
|
|
|
|
611
|
|
|
|
|
|
|
$svdvals = svdindexccs($u,$s,$v, $ptr,$rowids); |
612
|
|
|
|
|
|
|
$err_nz = ($nzvals-$svdvals)->pow(2)->sumover; |
613
|
|
|
|
|
|
|
|
614
|
|
|
|
|
|
|
|
615
|
|
|
|
|
|
|
|
616
|
|
|
|
|
|
|
=for bad |
617
|
|
|
|
|
|
|
|
618
|
|
|
|
|
|
|
svderror does not process bad values. |
619
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
620
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
|
622
|
|
|
|
|
|
|
=cut |
623
|
|
|
|
|
|
|
|
624
|
|
|
|
|
|
|
|
625
|
|
|
|
|
|
|
|
626
|
|
|
|
|
|
|
|
627
|
|
|
|
|
|
|
|
628
|
|
|
|
|
|
|
|
629
|
|
|
|
|
|
|
*svderror = \&PDL::svderror; |
630
|
|
|
|
|
|
|
|
631
|
|
|
|
|
|
|
|
632
|
|
|
|
|
|
|
|
633
|
|
|
|
|
|
|
|
634
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
635
|
|
|
|
|
|
|
=pod |
636
|
|
|
|
|
|
|
|
637
|
|
|
|
|
|
|
=head1 ACKNOWLEDGEMENTS |
638
|
|
|
|
|
|
|
|
639
|
|
|
|
|
|
|
Perl by Larry Wall. |
640
|
|
|
|
|
|
|
|
641
|
|
|
|
|
|
|
PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others. |
642
|
|
|
|
|
|
|
|
643
|
|
|
|
|
|
|
SVDLIBC by Dough Rohde. |
644
|
|
|
|
|
|
|
|
645
|
|
|
|
|
|
|
SVDPACKC by Michael Berry, Theresa Do, Gavin O'Brien, Vijay Krishna and Sowmini Varadhan. |
646
|
|
|
|
|
|
|
|
647
|
|
|
|
|
|
|
=cut |
648
|
|
|
|
|
|
|
|
649
|
|
|
|
|
|
|
##---------------------------------------------------------------------- |
650
|
|
|
|
|
|
|
=pod |
651
|
|
|
|
|
|
|
|
652
|
|
|
|
|
|
|
=head1 KNOWN BUGS |
653
|
|
|
|
|
|
|
|
654
|
|
|
|
|
|
|
Globals still lurk in the depths of SVDLIBC. |
655
|
|
|
|
|
|
|
|
656
|
|
|
|
|
|
|
=cut |
657
|
|
|
|
|
|
|
|
658
|
|
|
|
|
|
|
|
659
|
|
|
|
|
|
|
##--------------------------------------------------------------------- |
660
|
|
|
|
|
|
|
=pod |
661
|
|
|
|
|
|
|
|
662
|
|
|
|
|
|
|
=head1 AUTHOR |
663
|
|
|
|
|
|
|
|
664
|
|
|
|
|
|
|
Bryan Jurish Emoocow@cpan.orgE |
665
|
|
|
|
|
|
|
|
666
|
|
|
|
|
|
|
=head1 COPYRIGHT AND LICENSE |
667
|
|
|
|
|
|
|
|
668
|
|
|
|
|
|
|
Copyright (c) 2005-2015, Bryan Jurish. All rights reserved. |
669
|
|
|
|
|
|
|
|
670
|
|
|
|
|
|
|
This package is free software, and entirely without warranty. |
671
|
|
|
|
|
|
|
You may redistribute it and/or modify it under the same terms |
672
|
|
|
|
|
|
|
as Perl itself, either version 5.20.2 or any newer version of Perl 5 |
673
|
|
|
|
|
|
|
you have available. |
674
|
|
|
|
|
|
|
|
675
|
|
|
|
|
|
|
The SVDLIBC sources included in this distribution are themselves |
676
|
|
|
|
|
|
|
released under a BSD-like license. See the file |
677
|
|
|
|
|
|
|
F in the PDL-SVDLIBC source distribution |
678
|
|
|
|
|
|
|
for details. |
679
|
|
|
|
|
|
|
|
680
|
|
|
|
|
|
|
=head1 SEE ALSO |
681
|
|
|
|
|
|
|
|
682
|
|
|
|
|
|
|
perl(1), PDL(3perl), PDL::CCS(3perl), SVDLIBC documentation. |
683
|
|
|
|
|
|
|
|
684
|
|
|
|
|
|
|
=cut |
685
|
|
|
|
|
|
|
|
686
|
|
|
|
|
|
|
|
687
|
|
|
|
|
|
|
|
688
|
|
|
|
|
|
|
; |
689
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
|
691
|
|
|
|
|
|
|
|
692
|
|
|
|
|
|
|
# Exit with OK status |
693
|
|
|
|
|
|
|
|
694
|
|
|
|
|
|
|
1; |
695
|
|
|
|
|
|
|
|
696
|
|
|
|
|
|
|
|