- Join over
**1.5M+ people** - Join over
**100K+ communities** - Free
**without limits** - Create
**your own community**

I think they have something which implements it much much faster. I am trying to convert the code you have to

`C`

and still MATLAB is faster.

How do you compare the methods? The drop tolerance can have a large influence on performance, but MATLAB is using a locally weighted version, so that might make a difference.

I get lots of segmentation faults with Numba 0.50.1, I guess they broke something again. For now, downgrading the numba version should fix this, but I'll look into it.

I get

`TypeError: No matching definition for argument type(s) int64, array(int32, 1d, C), array(int64, 1d, C), array(int64, 1d, C), array(float64, 1d, C), array(int64, 1d, C), array(int64, 1d, C), float64, float64, int64`

on Numba 0.50.1. Any idea?

The first array should be of type `float64`

instead of `int32`

.

The first array should be of type float64 instead of int32.

I took your code as is. Or do you mean the input I send?

The first array should be of type float64 instead of int32.

I took your code as is. Or do you mean the input I send?

Your error message says that the second parameter to `_ichol`

is an array of type `int32`

, but it must be of type `float64`

instead, i.e. computing an incomplete Cholesky decomposition of an integer matrix is not supported.

@99991 , I will recheck my code.

Anyhow, have a look at: https://github.com/RoyiAvital/IncompleteCholeskyDecompositionThreshold for the*Run Time* comparison.

Anyhow, have a look at: https://github.com/RoyiAvital/IncompleteCholeskyDecompositionThreshold for the

I had a quick look at your code and found the following issues:

- MATLAB's implementation uses
`norm(A(j:end,j),1)*droptol`

instead of an absolute threshold. This might make a difference for the matrices you are testing with. You can compute the scaling factor by summing over`abs(Lij)`

here. Then check the number of non-zero elements of the results of both ichol implementations to see if it works correctly. - You are using quick sort, but the implementation from the standard library is not very fast because you have to pass a function pointer to it, so every comparison has the overhead of an indirect call. You can either use
`std::sort`

from C++, which is templated, so the code will be inlined. Alternatively, you can implement your own sorting method. The number of elements (at least for alpha matting) is usually small. In this case, sorting networks work well. For more elements, you can first use e.g. merge sort and then switch to sorting networks or insertion/selection sort for few elements. - You check if
`vD[jj]`

is less than zero, but ichol will also fail if the value is exactly zero since all column elements will be divided by it later. The check should be`<=`

instead of`<`

. - When
`discardThr == 0.0`

, the exact sparse Cholesky decomposition will be computed, which is a great way to check correctness, but your check here disallows that value.

Hi @99991 , Thanks for the pointers. At the begining I indeed used the

I am aware that MATLAB's thresholding approach depends on the row. But still, it is amazing how much faster MATLAB is.

Maybe there are better approaches for this decomposition? I saw this https://cerfacs.fr/wp-content/uploads/2016/03/scott.pdf.

`discardThr`

in the manner you mentioned to validate the code.I am aware that MATLAB's thresholding approach depends on the row. But still, it is amazing how much faster MATLAB is.

Maybe there are better approaches for this decomposition? I saw this https://cerfacs.fr/wp-content/uploads/2016/03/scott.pdf.

MATLAB is still faster (For sorting). Probably the overhead.

Hi there, thanks for making this library - it looks great! I've just been reading through GUCH20 [Multilevel Approach] and I'm wondering how I can measure foreground extraction error in the same way as the paper (with the adapted Rhemann et. al's metrics). I'm quite new to this field so any help would be greatly appreciated!

MATLAB is doing something very efficient. The

`WLS`

matrix should be very similar to your case and MATLAB's ICT is doing amazingly well there.
hey all! im using the simple pymatting example to extract the foreground of some images. the results are generally pretty good, but there is this odd halo region around the entire object that looks like artifacts of some sort. does anyone know what this could be from?

the image is a .png, quite high resolution (6000x4000), so im scaling the image down by a quarter when i load it with pymatting's

`load_image`

utility...i thought maybe the artifacts are from the resizing operation? but i tried other resolutions and the issue seems to persist
@MichaelWalczyk_twitter Impossible to tell without input image and trimap.

@ismailsalim The code can be easily adapted from here https://github.com/poppinace/indexnet_matting/tree/master/evaluation_code

Could it be that MATLAB uses AMD (Approximate Minimum Degree Permutation) - https://www.mathworks.com/help/matlab/ref/amd.html before applying the

Probably not as it means it will need to propagate the permutation to other functions as well.

`ichol()`

?Probably not as it means it will need to propagate the permutation to other functions as well.

@99991 When using

How exactly do I go about changing those parameters?

`estimate_alpha_cf()`

I get the following warning:```
PERFORMANCE WARNING:
Thresholded incomplete Cholesky decomposition failed due to insufficient positive-definiteness of matrix A with parameters:
discard_threshold = 1.000000e-04
shift = 0.000000e+00
Try decreasing discard_threshold or start with a larger shift
```

How exactly do I go about changing those parameters?

@99991 Also, am I right in thinking that it's not possible to use

`estimate_alpha_cf()`

with sparse scribbles as opposed to a dense trimap?
@ismailsalim You can either pass your own

`preconditioner`

function into the `estimate_alpha_cf`

function or you can use the expert level interface where you can call `ichol`

directly and pass different `discard_threshold`

and `shifts`

parameters. However, if you are using very sparse trimaps, the closed form alpha matting matrix is particularly hard to invert, in which case `ichol`

will have more trouble to succeed. Usually, sparse trimaps can be extended to denser trimaps during construction, for example by propagating the known and unknown pixels into regions with very similar color, like a paint bucket tool. Another option is to only trace rough outlines of objects and fill the interior and exterior. This way, the number of unknown pixels can be greatly reduced without too much effort. Or if you do not want to deal with stability issues and have less strict performance requirements, you can also use the `smoothed_aggregation`

preconditioner from pyamg as described here pymatting/pymatting#19
@RoyiAvital I had access to a MATLAB machine and their

`ichol`

implementation seemed to be slightly faster than our numba implementation, but then I translated the Python code to C and the performance difference went away. I believe that most of the performance difference is due to the choice of the threshold parameter, but I will try a WLS matrix for image colorization when I can find some time since I believe that in your benchmarks the performance difference was greatest for this kind of linear system. The performance also seemed bad with random matrices, but those matrices also have high fill-in, so a dense cholesky decomposition might be more appropriate in this case and so far I never came across such matrices in real world applications anyway.
@RoyiAvital Here is some extremely messy C code using the WLS matrix for image enhancement: https://github.com/99991/testing Unfortunately I will not have time to clean it up for several months due to teaching obligations. However, the performance of matlab and this code is already extremely close without too many optimizations, so maybe it is useful to you anyway.

@RoyiAvital These timings only include the time for

`ichol`

, no image loading or saving, so for MATLAB, it is the time from line 32 to 36.
I will try to wrap it as a MEX and compare it to see.

@99991 , I understand that

`relative_threshold`

(From https://github.com/99991/testing/blob/49de408bc9807e7b9090aba9318f4ce3d03eb2bb/ichol.c#L79) is to imitate MATLAB's tolerance parameter.
@99991 , I wrapped your code into MEX (See https://github.com/RoyiAvital/IncompleteCholeskyDecompositionThreshold, as

`RoyiIcholMex.c`

). It seems there is a bug with the relative threshold.