From octave-sources-request at bevo dot che dot wisc dot edu Thu Jan 11 09:24:15 2001 Subject: dot.m (compromise revision) From: Rolf Fabian To: "'octave-sources UWISC'" Date: Thu, 11 Jan 2001 16:20:46 +0100 Hi, Hope this upload is converging against a last revision concerning this topic I've got a reply by Dirk Laurie (see below) which sheds more llight on the problem. He observed a tendency, that more ,pure' orientated mathematical textbooks use the dot(x,y) definition which conjugates the SECOND argument, whereas more ,applied numerical' orientated linear algebra textbooks prefer complex conjugation of the FIRST argument (see below). Now, because Octave is more a numerical application plattform than ,pur', it's more natural' to give the default conjugation preference to the first argument ( which is also fully compatible to Matlab ). In other words, Dirk convinced me to change the default mode preference with respect to earlier revisions. (Best to delete them). Thanks Dirk! Furthermore I've obtained several valuable suggestions by Paul Kienzle concerning execution speedup, especially some hints how to work around the slow ,exist()' function. I've implemented them. Thanks Paul! ----------------- REVISED FUNCTION CHECK -------------- :) x =[ 1 + 0i 1 - 1i 1 + 1i ] :) y =[ 3 - 1i 3 + 1i 0 + 2i ] #by default we're in ,applied mathematican's world :)dot(x,y), dot(x.',y), dot(x,y.'), dot(x.',y.') #check all 4 orientation combinations ans = 7 + 5i #Ok this is same as Matlab output! ans = 7 + 5i ans = 7 + 5i ans = 7 + 5i :)global dot_arg_conj_prefer=2; #switch to ,pure' mathematician's world! :)dot(x,y), dot(x.',y), dot(x,y.'), dot(x.',y.') #check all 4 orientation combinations ans = 7 - 5i ans = 7 - 5i ans = 7 - 5i ans = 7 - 5i :)dot(x,y,3) #test Matlab 3rd arg comatibility! error: dot: multidimensional arrays not available for Octave ----------------- REVISED FUNCTION CHECK -------------- ---------------------------------------------------- Reply by Dirk Laurie: Dear Rolf: Sorry, I misread your posting and replied too hastily. Please accept my apologies for suggesting that you could do something so illogical. I personally have never used a textbook that conjugates the second argument, which is the reason why your point did not penetrate my thick skull. However, when I consulted my shelves, I found [3] Sheldon Axler, "Linear Algebra Done Right" Springer On the other hand, many standard textbooks conjugate the first argument. [4] Gilbert Strang, "Introduction to Linear Algebra" Wellesley-Cambridge [5] Golub and Van Loan, "Matrix Computations" Johns Hopkins [6] G.W. Stewart, "Matrix Algorithms" SIAM [7] J.H. Wilkinson, "The Algebraic Eigenvalue problem" Oxford I seem to detect a pattern in the difference. My reference [3], as well as the two you give > [1] S. Lipschutz, ,Linear Algebra' Schaum (English) > [2] Kowalski, ,Linear Algebra' DeGruyter (German) are books by "pure" mathematicians, and the other four by "numerical" mathematicians. Especially [7], first published in 1965 and often reprinted (my copy is a sixth printing dated 1978), has been very influential. Here is a theory that explains the difference. A pure mathematician introduces real dot products in a form like (A,B) = a_1 b_1 + a_2 b_2 + ... + a_n b_n which does not distinguish between A and B. Later, when complex numbers are introduced, it does not matter which you choose to conjugate as long as you stick to your choice. A numerical mathematician introduces matrix multiplication at an early stage, and identifies a vector with an nx1 matrix. Therefore the inner product is introduced as (A,B) = A^T B Now one is in any case going to generalize "transpose" of matrices to "Hermitian transpose", and it is natural to prefer the definition (A,B) = A^H B instead of A^T conj(B). In a setting where one identifies a vector with an 1xn matrix, of course, conjugating the second argument is the natural thing to do. This does however seem to be less usual in numerical linear algebra practice, Sorry, again, for impugning your good intentions. Dirk ---------------------------------------------------- -------------- REVISED FUNCTION SCRIPT ---------- ## -*- texinfo -*- ## at deftypefn {Function File} {} dot (@var{x}, @var{y}, @var{dim}) ## DOT PRODUCT of two vectors ## ##The dot (inner) product may be legally defined in two flavours: ##(A) at conj(var{x}(1))*@var{y}(1) + .. + @conj(var{x}(N))*@var{y}(N) ## This definition prefers complex conjugation of the first argument. ## It is mostly used by 'applied' linear algebra textbooks (and Matlab). ## Choosen to be the default method in the current function. ## ##(B) at var{x}(1)*conj(@var{y}(1)) + .. + @var{x}(N)*conj(@var{y}(N)) ## This alternate definition applies complex conjugation to the second ## argument and is preferred by more 'pure' mathematical textbooks. ## ##Both definitions are identical for vectors at var{x}, @var{y} in ##euclidian (real) space, but differ in unitary (complex) space. ##Switching between 'applied' and 'pure' mathematician's world ##may be controlled by setting-up a variable ## 'global dot_arg_conj_prefer' =1 (A) or =2 (B). ## ##The third argument at var{dim} is implemented for Matlab compatibility. ##It's usage invokes an error message that multidimensional arrays are ##not available for Octave. ## ##Simple properties ##dot(C* at var(x),@var(y)) == conj( C ) * dot(@var(x),@var(y)) (A) ## == C * dot( at var(x),@var(y)) (B) ##dot( at var(x),C*@var(y)) == C * dot(@var(x),@var(y)) (A) ## == conj( C ) * dot( at var(x),@var(y)) (B) ##dot( ) == norm( at var(x) )^2 (A ) == norm( @var(x) )^2 (A|B) ##dot( ) == conj( dot( at var(y),@var(x) ) ) (A ) == conj( dot( @var(y),@var(x) ) ) (A|B) ## ##The input vectors may pocess any combination of row/column ##orientation, but their length N must be the same. ## at end deftypefn ##AUTHOR Rolf Fabian Jan. 11, 2001 (rev.2) ## published under current GNU GENERAL PUBLIC LICENSE function z = dot(x,y,dim) if (nargin != 2) if nargin==3 error("dot: multidimensional arrays not available for Octave"); else usage ("dot( x,y,dim )"); endif endif global dot_arg_conj_prefer=1; DACP=dot_arg_conj_prefer; save_do_fortran_indexing = do_fortran_indexing; unwind_protect do_fortran_indexing = 1; if ( any(size(x)==1) && any(size(y)==1) && length(x)==length(y) ) if DACP==1 z = x(:)' * y(:); #conj 1st arg elseif DACP==2 z = x(:).' * conj( y(:) ); #conj 2nd arg endif else error ("dot: both arguments must be vectors of the same length"); endif unwind_protect_cleanup do_fortran_indexing = save_do_fortran_indexing; end_unwind_protect endfunction ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------