patchGNU Octave - Patches: patch #10067, [octave forge] (linear-algebra)...

 
 

patch #10067: [octave forge] (linear-algebra) Improvement of funm

Submitter:  Simon <sim79>
Submitted:  Fri 14 May 2021 10:23:10 PM UTC
   
 
Category:  Forge : other Priority:  3 - Low
Status:  None Privacy:  Public
Assigned to:  None Open/Closed:  Open
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 17 Jul 2022 05:31:06 PM UTC, comment #12: 

I attached the patch. The complete package was only intended for testing.
I agree with you that we should keep the current implementation of funm. So the patch adds a new function funm_sp and keeps funm. This also ensures backward compatibility.

(file #53448)

Simon <sim79>
Sat 16 Jul 2022 08:51:44 PM UTC, comment #11: 

Further to the funm.m discussion, you've uploaded a complete linear-algebra package.
But we rather have a patch or even better, a Mercurial changeset, wherein the new files are added and old files are dropped and that contains an explanation of exactly what things were changed. We cannot replace a complete package in our on-line repo, we can only handle changesets.

Also I wonder if that old function could be kept as 'ofunm.m' or so. My motive would be that it's quite a bit faster and requires less input. The trade-off is of course that it'll only work reliably for matrices with guaranteed distinct eigenvalues; in my use case (that I wrote funm.m for 22 years ago) that is invariably the case.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 06 Jul 2022 08:18:47 PM UTC, comment #10: 

Thanks for the explanation.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 06 Jul 2022 06:57:01 PM UTC, comment #9: 

Assume that input matrix A has the Schur decomposition A = QTQ* with upper triangular matrix T. The algorithm groups the eigenvalues of T into clusters. Eigenvalues in each cluster have similar absolute values. In order to create these clusters eigenvalues must be moved on the diagonal of T. E.g., eigenvalue T(i,i) has to be moved to an arbitrary position T(j,j).
In the original code on the homepage of Prof. Higham (Funm_files.tar) this moving is done with trexc, see files swap.c and swapping.m in the code. And in this paper (p. 475) Higham also suggests to use trexc.
ordschur moves a list of eingenvalues to the upper left of the diagonal of T. Quote from help string: 

> Reorders the real Schur factorization (U,S) obtained with the 'schur' function, so that selected eigenvalues appear in the upper left diagonal blocks of the quasi triangular Schur matrix.

But the algorithm does not always move eingenvalues to the upper left of the diagonal.

Simon <sim79>
Mon 04 Jul 2022 06:23:43 PM UTC, comment #8: 


> 3) ordschur is replaced by LAPACK-Routine xTREXC (see trexc.cc).

Just FTR ordschur has been implemented in Octave since 2015. In patch #8307, comment #6 it is mentioned that it is based on Lapack's 'trsen' rather than 'xTREXC'. Would you have any comment on that? I'm just asking, for the sake of easy maintainability, to be sure that calls to different routines that merely do the same thing are avoided.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 24 Dec 2021 10:06:00 PM UTC, comment #7: 

I included file funm.m from the RickOne-Repo into the linear-algebra package. Some minor changes to funm.m were necessary. List of changes:

1) funm is renamed to schur_parlett_funm since the RickOne-Repo and the linear-algebra package use different FUN parameters. In the RickOne-Repo parameter FUN has to be a function handle with two parameters x and k (for exp, sin, cos, log, sinh and cosh exceptions are possible). In the linear-algebra package the parameter FUN is a string or a function handle with only one parameter x. 
2) Help text of funm is converted into an octave docstring.
3) ordschur is replaced by LAPACK-Routine xTREXC (see trexc.cc).
4) In function swapping the line 

ind = cat (2,ind,{mdone:mdone+h(i)-1})

is replaced by

ind{i} = mdone:mdone+h (i)-1

since the original line causes a crash (zero indices are stored in cell array ind). The crash can be reproduced with

A = [ 4, 2, 0; 1, 4, 1; 1, 1, 4];
C = funm(A,@exp);

5) Additional BISTs.
6) Additional references in docstring.

(file #52564)

Simon <sim79>
Sat 22 May 2021 01:12:40 PM UTC, comment #6: 


> I'm more interested in a project in which ...

Fair enough.

> Did you write such a function handle for the Bessel function?

I can't remember (I probably did, but I'd have to dig it up). I used the modified Bessel function K for which there are fairly simple recurrence relations between derivatives.

ATM I have no time so I won't dig up the funm.m discussion. As to the GSOC project, I suppose it'll stall for another year or couple of years.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 21 May 2021 07:24:03 PM UTC, comment #5: 

I'm more interested in a project in which I can deal with programming and algorithms. I think this is not the case here since most of the work has been done.

Mudit Sharma is still active on the maintainers mailing list. In February 2021 he offered help related to the RickOne16 repo, see this message https://octave.1599824.n4.nabble.com/Re-GSoC-2021-application-td4698569.html. Maybe he can tell us what needs to be done before funm is included into the LA-package. And maybe he (or an other developer from the GSOC 16 thread) can assess whether an implementation in C++ or vectorization make sense.

I checked out the RickOne16 repo and tested funm with Octave 7.0.0. All tests were passed. I took a closer look at the funm implementation. By default it supports exp, sin, cos, cosh, sinh, and log. For log a specialized function logm_isst is called. For all other functions the caller has to pass a function handle that allows the evaluation of the k-th derivative.

@Philip: Did you write such a function handle for the Bessel function?  

The implementation makes extensive use of code from Funm_files.tar mentioned in comment #3. For example compare funm_atom.

Simon <sim79>
Tue 18 May 2021 07:19:06 PM UTC, comment #4: 

Yeah, prof. Higham's group has probably been subsidized my TMW so he can't share that code.

The code at https://github.com/RickOne16/matrix was made under supervision of another linear algebra expert, but wasn't fully finished; I don't know why. I think it's just a few things and then it works.

If you're interested you can help to get it into the linear-algebra package and replace the funm.m kludge I wrote more than 20 years ago (actually my very first contribution to Octave). I have other priorities (like adding functions to the packages I maintain, e.g., patch #10011 which has stalled for months now due to lack of time) but I can assist and guide you. IIRC it was quite easy to get it working.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 18 May 2021 06:32:04 PM UTC, comment #3: 

thanks for your reply. On the General Discussion thread is a link to Funm_files.tar (http://www.maths.manchester.ac.uk/~higham/NAMF/). I looked at this code but in the thread you also find the following statement by N.J. Higham:

> [...] I don't have anything distributable due to funm being copyrighted by The MathWorks. However, I think it's true to say that the state of the art Schur-Parlett algorithm is fully described in the union of my various papers - and in particular in my book [...]

I wasn't sure whether this statement includes Funm_files.tar. So I tried to write the code from scratch based on the pseudocode in the book. But there is still some resemblance to Funm_files.tar.

The code at https://github.com/RickOne16/matrix is certainly better than my implementation and should be used for future development. My approach was more "experimental" than "professional". In particular, I didn't think about vectorization, JIT etc..

Perhaps we can close this thread.

Simon <sim79>
Sun 16 May 2021 03:29:20 PM UTC, comment #2: 

I've just update the wiki to add those references.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 16 May 2021 03:23:04 PM UTC, comment #1: 

(Why was this thread "private" and can it only be seen when being logged in?
I've reset Privacy flag -> Public)

Thank you for your contribution.
The funm stanza on the Projects page you referred to unfortunately turns out to be a bit terse; some time ago the wiki has been reorganized and since then I lost track of this particular subject.
Quite some work has already been done, but that is hard to find if there's no mention of it and no indication of where to search. Apologies for that.
Here's another thread:
https://octave.1599824.n4.nabble.com/Re-GSOC-16-Improvements-to-sqrtm-logm-and-funm-td4675180.html
and a working funm.m (in .m format, not compiled code) is available here:
https://github.com/RickOne16/matrix

I've been using that code quite some time ago and it worked fine for my needs. TBH I can't remember if it worked OOTB or if I had to fix andor add and/or compile a few things here and there.
There seem(ed) to be plans to rewrite it all in C++ for performance but that has never materialized.

Compared to your version the other code accepts about any function in Octave's workspace which makes it a more generally usable function; for example, I used it with modified Bessel functions. But I'm no expert in matrix functions so I can't judge very well which version has the most "potential".

Philip Nienhuis <philipnienhuis>
Group Member
Fri 14 May 2021 10:23:10 PM UTC, original submission:  

Hello,

this (experimental) package contains the implementation of function funm with the Schur Parlett algorithm from N.J. Higham, Functions of Matrices, 2008. In order to avoid conflicts with the existing funm implementation the function is named funm2. There is a proposal on the Octave Projects page which is related to this package. Only exp,sin,cos,sinh and cosh are supported by default because their k-th derivatives can be calculated easily.

I tested the package with Octave 7.0.0 on Ubuntu 20.4.

Future TODO list
*) implement k-th derivative for functions other than exp,sin,cos,sinh and cosh
*) move declarations of LAPACK procedures in trexc.cc to file lo-lapack-proto.h
*) check whether funm2 can be used for functions whose Taylor series doesn't have
   an infinite radius of convergence.
*) move funm2 to linear algebra package

Simon

Simon <sim79>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #53448:  funm_sp_rev1.patch added by sim79 (28KiB - text/x-patch)
file #52564:  linear-algebra-2.2.4.tar.gz added by sim79 (398KiB - application/gzip)
file #51434:  funm2-1.0.0.tar.gz added by sim79 (19KiB - application/gzip)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by philipnienhuis
  • -email is unavailable- added by siko1056 (Updated the item)
  • -email is unavailable- added by sim79 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-07-17 sim79 Attached File- Added funm_sp_rev1.patch, #53448
    2021-12-24 sim79 Attached File- Added linear-algebra-2.2.4.tar.gz, #52564
    2021-05-17 philipnienhuis Carbon-Copy- Added philipnienhuis
    2021-05-16 philipnienhuis PrivacyPrivate Public
    2021-05-15 siko1056 SummaryImprovement of funm (linear algebra package) [octave forge] (linear-algebra) Improvement of funm
    2021-05-14 sim79 Attached File- Added funm2-1.0.0.tar.gz, #51434

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code