symamd is not comparable with OCTAVE/MATLAB

[VikasChidananda]

I have been trying to convert some code from MATLAB to Julia. But I observe that symamd operation is not very comparable to MATLAB's.
I am pasting the whole code snippet for comparison:

Julia Code

H   =   2;                                 #Height of the container
L   =   5;                                 #Length of the container

nx    =   10;
ny    =   7; 

dx  =   L/(nx-1);                         #Width of space step(x)
dy  =   H/(ny-1);                         #Width of space step(y)
dt   = 0.01

e       =   ones(nx-2);
i       =   ones(ny-2);

Tx      =   spdiagm(-1 => e[1:end-1], 0 => -2*e, 1 => e[1:end-1]);      
Ty      =   spdiagm(-1 => i[1:end-1], 0 => -2*i, 1 => i[1:end-1]);

Tx[1,1] =   -1; 
Tx[end,end]=-1;

Tt      =   kron(Ty/dy^2,   sparse(I, nx-2, nx-2))    +   kron(sparse(I, ny-2, ny-2),   Tx/dx^2);
N       =   (nx-2)*(ny-2)
Tt      =   sparse(I, N, N)     -   dt*Tt/Pe;
pt      =   symamd(Tt);

Output in Julia:

pt = [1, 8, 33, 40, 39, 32, 34, 25, 16, 7, 9, 2, 4, 36, 37, 5, 10, 15, 26, 31, 19, 21, 24, 23, 30, 14, 22, 29, 38, 13, 6, 28, 12, 20, 27, 35, 18, 17, 11, 3]

Output in OCTAVE:

debug> pt
pt =

 Columns 1 through 11:

    1    9    2   10   33   34   25   26   19    4    5

 Columns 12 through 22:

   36   37   21    8   16    7   15   40   39   32   31

 Columns 23 through 33:

   24   23   30   14   22   29   38   13    6   28   12

 Columns 34 through 40:

   20   27   35   18   17   11    3

Any help/suggestion/insight is greatly appreciated,
Thanks!

[amontoison]

Hi @VikasChidananda!
symamd computes a fill-in reducing permutation but it's based on an heuristic. Do you always have the same permutation in Julia and in OCTAVE / MATLAB?
What is the function that you called in OCTAVE / MATLAB?

[VikasChidananda]

Hi @amontoison ,

yes! I just checked for 5 runs, it always gives the same permutation for the above code in OCTAVE and Julia (as in the code output above)
The equivalent function in OCTAVE is also called symamd

[amontoison]

In AMD.jl, we directly call the C routines of SuiteSparse.
For MATLAB / OCTAVE they probably use the *.m files of the this folder.
Maybe the default options are different.

Does the fill-in of the factors is different after a factorization?

[VikasChidananda]

Sorry for the late reply.
I realized the above is not the best representation of minimal working example. In fact yes, the fill are not comparable either! Consider the following code:

Julia Code

using LinearAlgebra
using SparseArrays
using AMD

Re = 1e2
H   =   2;                                 #Height of the container
L   =   5;                                 #Length of the container

nx    =   10;
ny    =   7; 

dx  =   L/(nx-1);                         #Width of space step(x)
dy  =   H/(ny-1);                         #Width of space step(y)
dt   = 0.01

e       =   ones(nx-2);
i       =   ones(ny-1);
Ax      =   spdiagm(-1 => e[1:end-1], 0 => -2*e, 1 => e[1:end-1]);
Ay      =   spdiagm(-1 => i[1:end-1], 0 => -2*i, 1 => i[1:end-1]);     
Ay[1,1] =   -3;
Ay[end,end]=-3;

A       =   kron(Ay/dy^2, sparse(I, nx-2, nx-2))   +   kron(sparse(I, ny-1, ny-1),    Ax/dx^2);
N       =   (nx-2)*(ny-1)
Dv      =   sparse(I, N, N)            -   dt*A/Re;

pv      =   symamd(Dv);
F       =   lu(Dv[pv, pv]);
Lv, Uv  =   F.L, F.U;

OCTAVE Code

Re = 1e2
H   =   2;                                 %Height of the container
L   =   5;                                 %Length of the container
nx    =   10;
ny    =   7; 
dx  =   L/(nx-1);                         %Width of space step(x)
dy  =   H/(ny-1);                         %Width of space step(y)
dt   = 0.01
e=ones(nx-2,1);
i=ones(ny-1,1);
Ax=spdiags(e*[1 -2 1],-1:1,nx-2,nx-2);
Ay=spdiags(i*[1 -2 1],-1:1,ny-1,ny-1);
Ay(1,1)=-3;
Ay(end,end)=-3;
A=kron(Ay/dy^2,speye(nx-2)) + kron(speye(ny-1),Ax/dx^2);
Dv=speye((nx-2)*(ny-1))-dt*A/Re;
pv=symamd(Dv);
[Lv Uv]=lu(Dv(pv,pv));

and if we check: diag(Lv, -1) for both cases the fill-ins are different.
Julia

julia> diag(Lv, -1)
47-element SparseVector{Float64, Int64} with 37 stored entries:
  [1 ]  =  -0.000322815
  [2 ]  =  -2.89751e-7
  [3 ]  =  -0.000323105
  [5 ]  =  -0.000896997
  [6 ]  =  -2.89824e-7
  [7 ]  =  -0.000896997
  [8 ]  =  -2.90272e-7
  [9 ]  =  -0.000323105
  [10]  =  -2.90178e-7
  [14]  =  -0.000322919
  [15]  =  -2.89657e-7
  [18]  =  -0.000322815
        ⋮
  [34]  =  -5.80357e-7
  [35]  =  -5.80544e-7
  [36]  =  -1.04431e-7
  [37]  =  -8.06051e-7
  [38]  =  -5.79576e-7
  [39]  =  -8.06051e-7
  [40]  =  -5.80357e-7
  [42]  =  -0.000323106
  [43]  =  -5.80357e-7
  [44]  =  -4.67796e-13
  [45]  =  -5.80096e-7
  [46]  =  -7.30028e-22
  [47]  =  -7.23026e-10

OCTAVE

>> diag(Lv, -1)
ans =

Compressed Column Sparse (rows = 47, cols = 1, nnz = 38 [81%])

  (1, 1) -> -8.9700e-04
  (2, 1) -> -2.8992e-07
  (3, 1) -> -2.6006e-10
  (4, 1) -> -2.9018e-07
  (7, 1) -> -3.2292e-04
  (8, 1) -> -2.8966e-07
  (9, 1) -> -3.2321e-04
  (10, 1) -> -3.2321e-04
  (11, 1) -> -2.9018e-07
  (14, 1) -> -3.2292e-04
  (15, 1) -> -2.8966e-07
  (18, 1) -> -8.9700e-04
  (19, 1) -> -2.8992e-07
  (20, 1) -> -2.6006e-10
  (21, 1) -> -2.9018e-07
  (22, 1) -> -8.9780e-04
  (23, 1) -> -2.9018e-07
  (26, 1) -> -3.2292e-04
  (27, 1) -> -2.8966e-07
  (28, 1) -> -3.2321e-04
  (30, 1) -> -3.2292e-04
  (31, 1) -> -8.9700e-04
  (32, 1) -> -2.8992e-07
  (33, 1) -> -8.9700e-04
  (34, 1) -> -5.8036e-07
  (35, 1) -> -3.2321e-04
  (36, 1) -> -2.8137e-10
  (37, 1) -> -5.8036e-07
  (38, 1) -> -5.8010e-07
  (39, 1) -> -9.3454e-13
  (40, 1) -> -5.8036e-07
  (41, 1) -> -5.8010e-07
  (42, 1) -> -4.6717e-13
  (43, 1) -> -5.8036e-07
  (44, 1) -> -1.0446e-07
  (45, 1) -> -1.2125e-13
  (46, 1) -> -5.8036e-07
  (47, 1) -> -8.9780e-04

[amontoison]

Hi @VikasChidananda, can you post the number of non-zeros of the factor Lv? The goal of the permutation is to reduce the fill-in of the factors (sparsity of the factors). It's normal that we have different values in Lv because we have different permutation in Julia / Octave.

[VikasChidananda]

I understand. Thanks for clarifying.
Here is the output for both:

julia> nnz(Lv)
241

>> nnz(Lv)
ans = 242

[amontoison]

Ok, so the two permutations give almost the same fill-in. It's what in want in practice.
If you want to compare the reorderings, nnz is the most relevant metric.

[dpo]

I would expect the permutations to be identical though. SYMAMD is deterministic as far as I know and only depends on the sparsity pattern. @VikasChidananda Have you confirmed that the latter is identical in Julia and
octave (e.g., by writing one of the matrices to file and reading it in at the other end)?

[dpo]

One possible difference is that we don't set default options explicitly:

• https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L90 vs
• https://github.com/JuliaSmoothOptimizers/AMD.jl/blob/main/src/COLAMD.jl#L169
  (should be the same but we should double check).

Another is these options: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L99

Also we don't return the stats array. It could help debug the issue.

finally, I notice that we don't call the "report" function in case of failure: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/dev/CCOLAMD/MATLAB/csymamdmex.c#L163

[VikasChidananda]

I would expect the permutations to be identical though. SYMAMD is deterministic as far as I know and only depends on the sparsity pattern. @VikasChidananda Have you confirmed that the latter is identical in Julia and octave (e.g., by writing one of the matrices to file and reading it in at the other end)?

Yes, I checked the sparsity pattern of Matrices from both. They are identical, with slight difference being OCTAVE stores the rounded values (until 5 decimal points) and Julia stores until (20 decimal points). The permutations of both matrices (let's call it octave_Dv, Julia_Dv) are the same in Julia. It differs from that of OCTAVE.

[dpo]

@VikasChidananda I wonder if #61 fixes your issue.