Mercurial > octave

comparison libinterp/corefcn/gsvd.cc @ 22235:63b41167ef1e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gsvd: new function imported from Octave-Forge linear-algebra package. * libinterp/corefcn/gsvd.cc: New function to the interpreter. Imported from the linear-algebra package. * CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, dbleGSVD.h: new classes imported from the linear-algebra package to compute gsvd of Matrix and ComplexMatrix. * liboctave/operators/mx-defs.h, liboctave/operators/mx-ext.h: use new classes. * libinterp/corefcn/module.mk, liboctave/numeric/module.mk: Add to the * scripts/help/__unimplemented__.m: Remove "gsvd" from list. * doc/interpreter/linalg.txi: Add to manual. build system. * NEWS: Add function to list of new functions for 4.2.
author Barbara Locsi <locsi.barbara@gmail.com>
date Thu, 04 Aug 2016 07:50:31 +0200
parents
children 065a44375723
comparison
equal deleted inserted replaced
22234:66dd260512a4 22235:63b41167ef1e
1 // Copyright (C) 1996, 1997 John W. Eaton
2 // Copyright (C) 2006, 2010 Pascal Dupuis <Pascal.Dupuis@uclouvain.be>
3 //
4 // This program is free software; you can redistribute it and/or modify it under
5 // the terms of the GNU General Public License as published by the Free Software
6 // Foundation; either version 3 of the License, or (at your option) any later
7 // version.
8 //
9 // This program is distributed in the hope that it will be useful, but WITHOUT
10 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 // details.
13 //
14 // You should have received a copy of the GNU General Public License along with
15 // this program; if not, see <http://www.gnu.org/licenses/>.
16
17 #ifdef HAVE_CONFIG_H
18 # include <config.h>
19 #endif
20
21 #include "defun.h"
22 #include "error.h"
23 #include "gripes.h"
24 #include "pr-output.h"
25 #include "utils.h"
26 #include "ovl.h"
27
28 #include "CmplxGSVD.h"
29 #include "dbleGSVD.h"
30
31 DEFUN (gsvd, args, nargout,
32 doc: /* -*- texinfo -*-
33 @deftypefn {Loadable Function} {@var{s} =} gsvd (@var{a}, @var{b})
34 @deftypefnx {Loadable Function} {[@var{u}, @var{v}, @var{c}, @var{s}, @var{x} [, @var{r}]] =} gsvd (@var{a}, @var{b})
35 @cindex generalised singular value decomposition
36 Compute the generalised singular value decomposition of (@var{a}, @var{b}):
37 @iftex
38 @tex
39 $$
40 U^H A X = [I 0; 0 C] [0 R]
41 V^H B X = [0 S; 0 0] [0 R]
42 C*C + S*S = eye(columns(A))
43 I and 0 are padding matrices of suitable size
44 R is upper triangular
45 $$
46 @end tex
47 @end iftex
48 @ifinfo
49
50 @example
51 u' * a * x = [I 0; 0 c] * [0 r]
52 v' * b * x = [0 s; 0 0] * [0 r]
53 c * c + s * s = eye(columns(a))
54 I and 0 are padding matrices of suitable size
55 r is upper triangular
56 @end example
57 @end ifinfo
58
59 The function @code{gsvd} normally returns the vector of generalised singular
60 values
61 @iftex
62 @tex
63 diag(C)./diag(S).
64 @end tex
65 @end iftex
66 @ifinfo
67 diag(r)./diag(s).
68 @end ifinfo
69 If asked for five return values, it computes
70 @iftex
71 @tex
72 $U$, $V$, and $X$.
73 @end tex
74 @end iftex
75 @ifinfo
76 U, V, and X.
77 @end ifinfo
78 With a sixth output argument, it also returns
79 @iftex
80 @tex
81 R,
82 @end tex
83 @end iftex
84 @ifinfo
85 r,
86 @end ifinfo
87 The common upper triangular right term. Other authors, like S. Van Huffel,
88 define this transformation as the simulatenous diagonalisation of the
89 input matrices, this can be achieved by multiplying
90 @iftex
91 @tex
92 X
93 @end tex
94 @end iftex
95 @ifinfo
96 x
97 @end ifinfo
98 by the inverse of
99 @iftex
100 @tex
101 [I 0; 0 R].
102 @end tex
103 @end iftex
104 @ifinfo
105 [I 0; 0 r].
106 @end ifinfo
107
108 For example,
109
110 @example
111 gsvd (hilb (3), [1 2 3; 3 2 1])
112 @end example
113
114 @noindent
115 returns
116
117 @example
118 ans =
119
120 0.1055705
121 0.0031759
122 @end example
123
124 @noindent
125 and
126
127 @example
128 [u, v, c, s, x, r] = gsvd (hilb (3), [1 2 3; 3 2 1])
129 @end example
130
131 @noindent
132 returns
133
134 @example
135 u =
136
137 -0.965609 0.240893 0.097825
138 -0.241402 -0.690927 -0.681429
139 -0.096561 -0.681609 0.725317
140
141 v =
142
143 -0.41974 0.90765
144 -0.90765 -0.41974
145
146 c =
147
148 0.10499 0.00000
149 0.00000 0.00318
150
151 s =
152 0.99447 0.00000
153 0.00000 0.99999
154 x =
155
156 0.408248 0.902199 0.139179
157 -0.816497 0.429063 -0.386314
158 0.408248 -0.044073 -0.911806
159
160 r =
161 -0.14093 -1.24345 0.43737
162 0.00000 -3.90043 2.57818
163 0.00000 0.00000 -2.52599
164
165 @end example
166
167 The code is a wrapper to the corresponding Lapack dggsvd and zggsvd routines.
168
169 @end deftypefn */)
170 {
171 octave_value_list retval;
172
173 int nargin = args.length ();
174
175 if (nargin < 2 || nargin > 2 || (nargout > 1 && (nargout < 5 || nargout > 6)))
176 {
177 print_usage ();
178 return retval;
179 }
180
181 octave_value argA = args(0), argB = args(1);
182
183 octave_idx_type nr = argA.rows ();
184 octave_idx_type nc = argA.columns ();
185
186 // octave_idx_type nn = argB.rows ();
187 octave_idx_type np = argB.columns ();
188
189 if (nr == 0 || nc == 0)
190 {
191 if (nargout == 5)
192 retval = ovl (identity_matrix (nc, nc), identity_matrix (nc, nc),
193 Matrix (nr, nc), identity_matrix (nr, nr),
194 identity_matrix (nr, nr));
195 else if (nargout == 6)
196 retval = ovl (identity_matrix (nc, nc), identity_matrix (nc, nc),
197 Matrix (nr, nc), identity_matrix (nr, nr),
198 identity_matrix (nr, nr),
199 identity_matrix (nr, nr));
200 else
201 retval = ovl (Matrix (0, 1));
202 }
203 else
204 {
205 if ((nc != np))
206 {
207 print_usage ();
208 return retval;
209 }
210
211 GSVD::type type = ((nargout == 0 || nargout == 1)
212 ? GSVD::sigma_only
213 : (nargout > 5) ? GSVD::std : GSVD::economy );
214
215 if (argA.is_real_type () && argB.is_real_type ())
216 {
217 Matrix tmpA = argA.matrix_value ();
218 Matrix tmpB = argB.matrix_value ();
219
220 if (! error_state)
221 {
222 if (tmpA.any_element_is_inf_or_nan ())
223 {
224 error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values");
225 return retval;
226 }
227
228 if (tmpB.any_element_is_inf_or_nan ())
229 {
230 error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values");
231 return retval;
232 }
233
234 GSVD result (tmpA, tmpB, type);
235
236 // DiagMatrix sigma = result.singular_values ();
237
238 if (nargout == 0 || nargout == 1)
239 {
240 DiagMatrix sigA = result.singular_values_A ();
241 DiagMatrix sigB = result.singular_values_B ();
242 for (int i = sigA.rows() - 1; i >=0; i--)
243 sigA.dgxelem(i) /= sigB.dgxelem(i);
244 retval = ovl (sigA.diag());
245 }
246 else
247 {
248 if (nargout > 5)
249 retval = ovl (result.left_singular_matrix_A (),
250 result.left_singular_matrix_B (),
251 result.singular_values_A (),
252 result.singular_values_B (),
253 result.right_singular_matrix (),
254 result.R_matrix ());
255 else
256 retval = ovl (result.left_singular_matrix_A (),
257 result.left_singular_matrix_B (),
258 result.singular_values_A (),
259 result.singular_values_B (),
260 result.right_singular_matrix ());
261 }
262 }
263 }
264 else if (argA.is_complex_type () || argB.is_complex_type ())
265 {
266 ComplexMatrix ctmpA = argA.complex_matrix_value ();
267 ComplexMatrix ctmpB = argB.complex_matrix_value ();
268
269 if (! error_state)
270 {
271 if (ctmpA.any_element_is_inf_or_nan ())
272 {
273 error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values");
274 return retval;
275 }
276 if (ctmpB.any_element_is_inf_or_nan ())
277 {
278 error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values");
279 return retval;
280 }
281
282 ComplexGSVD result (ctmpA, ctmpB, type);
283
284 // DiagMatrix sigma = result.singular_values ();
285
286 if (nargout == 0 || nargout == 1)
287 {
288 DiagMatrix sigA = result.singular_values_A ();
289 DiagMatrix sigB = result.singular_values_B ();
290 for (int i = sigA.rows() - 1; i >=0; i--)
291 sigA.dgxelem(i) /= sigB.dgxelem(i);
292 retval = ovl (sigA.diag());
293 }
294 else
295 {
296 if (nargout > 5)
297 retval = ovl (result.left_singular_matrix_A (),
298 result.left_singular_matrix_B (),
299 result.singular_values_A (),
300 result.singular_values_B (),
301 result.right_singular_matrix (),
302 result.R_matrix ());
303 else
304 retval = ovl (result.left_singular_matrix_A (),
305 result.left_singular_matrix_B (),
306 result.singular_values_A (),
307 result.singular_values_B (),
308 result.right_singular_matrix ());
309 }
310 }
311 }
312 else
313 {
314 gripe_wrong_type_arg ("gsvd", argA);
315 gripe_wrong_type_arg ("gsvd", argB);
316 return retval;
317 }
318 }
319
320 return retval;
321 }
322
323 /*
324 %# a few tests for gsvd.m
325 %!shared A, A0, B, B0, U, V, C, S, X, R, D1, D2
326
327 %! A0=randn(5, 3); B0=diag([1 2 4]);
328 %! A = A0; B = B0;
329 %! # disp('Full rank, 5x3 by 3x3 matrices');
330 %! # disp([rank(A) rank(B) rank([A' B'])]);
331 %! [U, V, C, S, X, R] = gsvd(A, B);
332 %! D1 = zeros(5, 3); D1(1:3, 1:3) = C;
333 %! D2 = S;
334 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
335 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
336 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
337
338 %! # disp('A 5x3 full rank, B 3x3 rank deficient');
339 %! B(2, 2) = 0;
340 %! [U, V, C, S, X, R] = gsvd(A, B);
341 %! D1 = zeros(5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C;
342 %! D2 = [zeros(2, 1) S; zeros(1, 3)];
343 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
344 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
345 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
346
347 %! # disp('A 5x3 rank deficient, B 3x3 full rank');
348 %! B = B0;
349 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
350 %! [U, V, C, S, X, R] = gsvd(A, B);
351 %! D1 = zeros(5, 3); D1(1:3, 1:3) = C;
352 %! D2 = S;
353 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
354 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
355 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
356
357 %! # disp("A 5x3, B 3x3, [A' B'] rank deficient");
358 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
359 %! [U, V, C, S, X, R] = gsvd(A, B);
360 %! D1 = zeros(5, 2); D1(1:2, 1:2) = C;
361 %! D2 = [S; zeros(1, 2)];
362 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
363 %!assert(norm((U'*A*X)-D1*[zeros(2, 1) R]) <= 1e-6)
364 %!assert(norm((V'*B*X)-D2*[zeros(2, 1) R]) <= 1e-6)
365
366 %! # now, A is 3x5
367 %! A = A0.'; B0=diag([1 2 4 8 16]); B = B0;
368 %! # disp('Full rank, 3x5 by 5x5 matrices');
369 %! # disp([rank(A) rank(B) rank([A' B'])]);
370
371 %! [U, V, C, S, X, R] = gsvd(A, B);
372 %! D1 = [C zeros(3,2)];
373 %! D2 = [S zeros(3,2); zeros(2, 3) eye(2)];
374 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
375 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
376 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
377
378 %! # disp('A 5x3 full rank, B 5x5 rank deficient');
379 %! B(2, 2) = 0;
380 %! # disp([rank(A) rank(B) rank([A' B'])]);
381
382 %! [U, V, C, S, X, R] = gsvd(A, B);
383 %! D1 = zeros(3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C;
384 %! D2 = zeros(5, 5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye(2);
385 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
386 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
387 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
388
389 %! # disp('A 3x5 rank deficient, B 5x5 full rank');
390 %! B = B0;
391 %! A(3, :) = 2*A(1, :) - A(2, :);
392 %! # disp([rank(A) rank(B) rank([A' B'])]);
393 %! [U, V, C, S, X, R] = gsvd(A, B);
394 %! D1 = zeros(3, 5); D1(1:3, 1:3) = C;
395 %! D2 = zeros(5, 5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye(2);
396 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
397 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
398 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
399
400 %! # disp("A 5x3, B 5x5, [A' B'] rank deficient");
401 %! A = A0.'; B = B0.';
402 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
403 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
404 %! # disp([rank(A) rank(B) rank([A' B'])]);
405 %! [U, V, C, S, X, R]=gsvd(A, B);
406 %! D1 = zeros(3, 4); D1(1:3, 1:3) = C;
407 %! D2 = eye(4); D2(1:3, 1:3) = S; D2(5,:) = 0;
408 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
409 %!assert(norm((U'*A*X)-D1*[zeros(4, 1) R]) <= 1e-6)
410 %!assert(norm((V'*B*X)-D2*[zeros(4, 1) R]) <= 1e-6)
411
412 %! A0 = A0 +j * randn(5, 3); B0 = B0=diag([1 2 4]) + j*diag([4 -2 -1]);
413 %! A = A0; B = B0;
414 %! # disp('Complex: Full rank, 5x3 by 3x3 matrices');
415 %! # disp([rank(A) rank(B) rank([A' B'])]);
416 %! [U, V, C, S, X, R] = gsvd(A, B);
417 %! D1 = zeros(5, 3); D1(1:3, 1:3) = C;
418 %! D2 = S;
419 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
420 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
421 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
422
423 %! # disp('Complex: A 5x3 full rank, B 3x3 rank deficient');
424 %! B(2, 2) = 0;
425 %! # disp([rank(A) rank(B) rank([A' B'])]);
426
427 %! [U, V, C, S, X, R] = gsvd(A, B);
428 %! D1 = zeros(5, 3); D1(1, 1) = 1; D1(2:3, 2:3) = C;
429 %! D2 = [zeros(2, 1) S; zeros(1, 3)];
430 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
431 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
432 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
433
434 %! # disp('Complex: A 5x3 rank deficient, B 3x3 full rank');
435 %! B = B0;
436 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
437 %! # disp([rank(A) rank(B) rank([A' B'])]);
438 %! [U, V, C, S, X, R] = gsvd(A, B);
439 %! D1 = zeros(5, 3); D1(1:3, 1:3) = C;
440 %! D2 = S;
441 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
442 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
443 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
444
445 %! # disp("Complex: A 5x3, B 3x3, [A' B'] rank deficient");
446 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
447 %! # disp([rank(A) rank(B) rank([A' B'])]);
448 %! [U, V, C, S, X, R] = gsvd(A, B);
449 %! D1 = zeros(5, 2); D1(1:2, 1:2) = C;
450 %! D2 = [S; zeros(1, 2)];
451 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
452 %!assert(norm((U'*A*X)-D1*[zeros(2, 1) R]) <= 1e-6)
453 %!assert(norm((V'*B*X)-D2*[zeros(2, 1) R]) <= 1e-6)
454
455 %! # now, A is 3x5
456 %! A = A0.'; B0=diag([1 2 4 8 16])+j*diag([-5 4 -3 2 -1]);
457 %! B = B0;
458 %! # disp('Complex: Full rank, 3x5 by 5x5 matrices');
459 %! # disp([rank(A) rank(B) rank([A' B'])]);
460 %! [U, V, C, S, X, R] = gsvd(A, B);
461 %! D1 = [C zeros(3,2)];
462 %! D2 = [S zeros(3,2); zeros(2, 3) eye(2)];
463 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
464 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
465 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
466
467 %! # disp('Complex: A 5x3 full rank, B 5x5 rank deficient');
468 %! B(2, 2) = 0;
469 %! # disp([rank(A) rank(B) rank([A' B'])]);
470 %! [U, V, C, S, X, R] = gsvd(A, B);
471 %! D1 = zeros(3, 5); D1(1, 1) = 1; D1(2:3, 2:3) = C;
472 %! D2 = zeros(5,5); D2(1:2, 2:3) = S; D2(3:4, 4:5) = eye(2);
473 %!assert(norm(diag(C).^2+diag(S).^2 - ones(2, 1)) <= 1e-6)
474 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
475 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
476
477 %! # disp('Complex: A 3x5 rank deficient, B 5x5 full rank');
478 %! B = B0;
479 %! A(3, :) = 2*A(1, :) - A(2, :);
480 %! # disp([rank(A) rank(B) rank([A' B'])]);
481 %! [U, V, C, S, X, R] = gsvd(A, B);
482 %! D1 = zeros(3, 5); D1(1:3, 1:3) = C;
483 %! D2 = zeros(5,5); D2(1:3, 1:3) = S; D2(4:5, 4:5) = eye(2);
484 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
485 %!assert(norm((U'*A*X)-D1*R) <= 1e-6)
486 %!assert(norm((V'*B*X)-D2*R) <= 1e-6)
487
488 %! # disp("Complex: A 5x3, B 5x5, [A' B'] rank deficient");
489 %! A = A0.'; B = B0.';
490 %! A(:, 3) = 2*A(:, 1) - A(:, 2);
491 %! B(:, 3) = 2*B(:, 1) - B(:, 2);
492 %! # disp([rank(A) rank(B) rank([A' B'])]);
493 %! [U, V, C, S, X, R]=gsvd(A, B);
494 %! D1 = zeros(3, 4); D1(1:3, 1:3) = C;
495 %! D2 = eye(4); D2(1:3, 1:3) = S; D2(5,:) = 0;
496 %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6)
497 %!assert(norm((U'*A*X)-D1*[zeros(4, 1) R]) <= 1e-6)
498 %!assert(norm((V'*B*X)-D2*[zeros(4, 1) R]) <= 1e-6)
499
500 */