# Extract of left and right vectors of multidimensional SVD

**URL:** https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844
**Category:** General Discussion
**Created:** [March 9, 2021, 2:16pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844 "2021-03-09T14:16:41Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 9, 2021, 2:16pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/1 "2021-03-09T14:16:42Z")

</div>

Hi,

I try to adapt the example “mf-2d-SLEPc.edp” to compute singular values for a multi-dimensional problem. Namely, I want to compute the singular values of the operator `MF`

> `Mat<complex> MF(Qu, Qf, prodFunc, transpose = prodFuncTranspose);`

where Qu is a matrix created on PuPuPp and Qf is a matrix created on PuPu where:

> `func PuPuPp=[P2, P2, P1];`  
> `func PuPu=[P2, P2];`

The associated spaces are:

> `fespace Uvvp(Th, PuPuPp);`  
> `fespace Uvv(Th, PuPu);`

I wrote the following lines to compute the singular values and extract the left and right vectors (the example helmholtz-2d-SLEPc-complex.edp helped me here):

> `int nsv = 4;`  
> `real[int] values(nsv);`  
> `Uvv<complex>[int] [rsolx,rsoly](nsv);`  
> `Uvvp<complex>[int] [lsolx, lsoly, lsolp](nsv);`  
> ` int nconv = SVDSolve(MF, sparams = "-svd_largest -svd_view_values -svd_type cyclic -svd_nsv "+nsv, values = values, lvectors = lsolx, rvectors = rsolx);`

SLEPc computes singular values that I can access in the array `values`. However the left and right vectors are all zero. I can observe that for instance by adding the following lines:

> for (int idSv = 0; idSv \< nconv; ++idSv)  
> {  
> cout \<\< "MAX " \<\< (rsolx[idSv]).linfty \<\< " / " \<\< (lsolx[idSv]).linfty \<\< endl;  
> }

Is there something wrong in the syntax to call `SVDSolve` to define `lvectors` and `rvectors`? Or is there something wrong about how I access them one the computation is done?

Best,

Lucas

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 9, 2021, 2:28pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/2 "2021-03-09T14:28:28Z")

</div>

Could you please provide a full [MWE](https://en.wikipedia.org/wiki/Minimal_working_example)? You can send it in private if it’s closed-source code. Thanks.

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 9, 2021, 2:39pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/3 "2021-03-09T14:39:03Z")

</div>

Actually by I modified the example “mf-2d-SLEPc.edp” (see attachment) [mf-2d-SLEPc.edp](https://community.freefem.org/uploads/short-url/vQjrjMVUAxeIdjjsgQ60NWc47me.edp) (2.3 KB) and I observe the same issue. So it is not due to the fact that my example is multi-dimensional. I guess I am doing something wrong when accessing the vectors once SVDSolve was called.

Thanks for your help

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 9, 2021, 2:50pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/4 "2021-03-09T14:50:03Z")

</div>

I see. There is no domain decomposition attached to `MF`, so you cannot use output vectors stored in FreeFEM numbering (`lvectors` and `rvectors`). Instead, you need to use vectors stored in PETSc numbering (`larray` and `rarray`), and if you need to go back to FreeFEM numbering, use `ChangeNumbering(..., inverse = true)`. Does that make sense to you?

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 9, 2021, 3:34pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/5 "2021-03-09T15:34:25Z")

</div>

Not completely!

> There is no domain decomposition attached to `MF`

MF is not indirectly aware of the decomposition made by `buildDmesh(Th);` ? Is there a way to attach the domain decomposition to `MF` ?

I roughly understood from “mf-2d-SLEPc.edp” how to use `ChangeNumbering`. I understand that I should replace `lvectors = lvec` by `larray = lvec` where `lvec` is now `real[int,int]` of size (Vh.ndof, 10) (as in [Composite eigenvectors with SLEPc](https://community.freefem.org/t/composite-eigenvectors-with-slepc/231))

Best,

Lucas

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 9, 2021, 4:13pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/6 "2021-03-09T16:13:30Z")

</div>

> MF is not indirectly aware of the decomposition

Correct. `C` is, not `MF`.

> I tried `real[int][int]` but it doest not seem to be correct.

What you need is `real[int, int`] (bidimensional array). `real[int][int]` is for arrays of arrays, but here, all columns (or rows for that matter) are of the same size.

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 9, 2021, 5:44pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/7 "2021-03-09T17:44:30Z")

</div>

> Correct. `C` is, not `MF` .

But then using `lvectors` on `C` should work? It does not to be the case as you can have observed from the above attached file.

I would like to extend the example “mf-2d-SLEPc.edp” to show how to use `larray`/`rarray` and then push it on the FreeFem Git as I think it adds some valuable explanation on how to use `SVDSolve`. I would like to add three asserts:

- same left vectors with and without matrix-free resolution
- same right vectors with and without matrix-free resolution
- matrix x (left vector) = (singular value) x (right vector)

Please find enclosed what I did so far: [mf-2d-SLEPc.edp](https://community.freefem.org/uploads/short-url/7CpVtNm4Zi78jdak9R2w2TyrZ2V.edp) (2.9 KB) . It is not functional yet: the comparison of left vectors fails. Moreover I struggle to change the numbering of the right vectors. If I understand correctly the matrix `B` should be used to change the numbering of the right vectors, but the re-numbering only works if I use the matrix `A`. Could you, please, help me to correct this?

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 9, 2021, 8:56pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/8 "2021-03-09T20:56:54Z")

</div>

Good catch! Should read `(*rarray)(':', i) = cpy;`, not `(*array)(':', i) = cpy;` in [SLEPc-code.hpp](https://github.com/FreeFem/FreeFem-sources/blob/f9e94258046aed5de8fc269e40747bea8e66400c/plugin/mpi/SLEPc-code.hpp#L464).

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 10, 2021, 10:26am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/9 "2021-03-10T10:26:01Z")

</div>

I compiled FreeFem with the proposed correction in SLEPc-code.hpp. The assert on the left vectors now successes. However, it still fail to handle the right vectors, [mf-2d-SLEPc.edp](https://community.freefem.org/uploads/short-url/zlNAhA45OtuyKAY9mcygCKegNgL.edp) (3.1 KB)  
I do not understand what I do wrong but I can not renumber the right vectors using the matrix `B`.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 10, 2021, 10:36am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/10 "2021-03-10T10:36:34Z")

</div>

Another good catch! Sorry, I’m not using `SVDSolve()` much, so thank you for your careful investigation. [This](https://github.com/FreeFem/FreeFem-sources/blob/c68fe0b/plugin/mpi/SLEPc-code.hpp#L399-L402) part of `SLEPc-code.hpp` is wrong, `rarray` can’t be of the same size as `array` in the rectangular case. That would explain why the `ChangeNumbering()` with `A` works (correct dimension) while it doesn’t with `B` (incorrect dimension). But you are correct, it should be the other way around.

I’m AFK until the early afternoon, but I’ll try to fix this later today. Thank you for bearing with me 😟.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 10, 2021, 3:56pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/11 "2021-03-10T15:56:16Z")

</div>

It should be good [now](https://github.com/FreeFem/FreeFem-sources/commit/d32eb3071c3cf3b56f82d930f0e36352368a7790#diff-af656507d13e13fb85236878286bb3ad84b78bb17a9a8fd41365bdeb7018b764). I’ve fixed [mf-2d-SLEPc.edp](https://community.freefem.org/uploads/short-url/tHQyOOdrpbWa8Y5jVZNuxVxw4UV.edp) (3.1 KB) accordingly, is it alright on your side if I commit this modified file on the FreeFEM official repository?

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 10, 2021, 4:24pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/12 "2021-03-10T16:24:23Z")

</div>

Thanks for the quick code correction!

Yes please commit the modified example 🙂

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 11, 2021, 8:20am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/13 "2021-03-11T08:20:57Z")

</div>

Actually there is one more question from my side. The last assert corresponds to check that C v = sigma u in the PETSc world. If I understand correctly, this should be true as well in the FreeFem world, isn’t it ?

Such test in FreeFem world would read:  
` temp1 = MF * vMF[];`  
`temp2 = valuesMF[idSV] * uMF[];`  
`temp2 -= temp1;`  
`assert(temp2.linfty < 1.0e-6);`

I tried that but that fails, did I do something wrong?

Best,

Lucas

**EDIT:**  
The following assert (that fails) convinced me that something is wrong when performing the matrix vector multiplication in the FreeFem world:

`MatMult(MF, rvecMF(:,idSV), temp2);`  
`ChangeNumbering(A, temp1, temp2, inverse = true, exchange = true);`  
`temp2 = MF * vMF[];`  
`temp2 -= temp1;`  
`assert(temp2.linfty < 1.e-6);`

It is not allowed to do `MF * vMF[]` ?

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 11, 2021, 8:46am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/14 "2021-03-11T08:46:17Z")

</div>

~~Can’t do `MF *` because `MF` is rectangular and I haven’t done the necessary plumbing to handle this case. Must use `MatMult` instead. If you had a FreeFEM compiled with debugging turned on, you would see an assert error when doing either `temp1 = MF * vMF[];` or `temp2 = MF * vMF[];`, whatever comes next is garbage.~~ (maybe not that simple, let me check and I’ll get back at you)

Edit: you forgot to resize `temp2` before `MF *`, the following seems to work.

```auto
MatMult(MF, rvecMF(:,idSV), temp2);
ChangeNumbering(A, temp1, temp2, inverse = true, exchange = true);
temp2.resize(uMF[].n); // should be done automatically, but necessary plumbing not done yet...
temp2 = MF * vMF[];
temp2 -= temp1;
assert(temp2.linfty < 1.e-6);

```

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 11, 2021, 9:31am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/15 "2021-03-11T09:31:42Z")

</div>

I see, thank you for the explanation 🙂

And just a remark for those like that would be interested (like me ^^), it is possible to stay in the “PETSc world” in prodFunc. In such case, prodFunc would read:

> func real[int] prodFunc(real[int]& up) {  
> real[int] u(A.n);  
> MatMultTranspose(RPETSc, up, u); // where RPETSc is defined by Mat RPETSc(B, A, R);  
> MatMult(A, u, up);  
> return up;  
> }

I guess this should improve the performance as we avoid two re-numbering.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 11, 2021, 9:58am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/16 "2021-03-11T09:58:18Z")

</div>

Indeed, it’s always best to stick as much as possible to PETSc numbering. This was merely to show the intent of `ChangeNumbering()` & friends 😄

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 16, 2021, 5:21pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/17 "2021-03-16T17:21:09Z")

</div>

Hi!

It seems that the correction you proposed broke something when using `EPSolve( ... , array = vec);`. I came to this conclusion by comparing the result of the same script executed with the v4.8 and the develop branch. In the last case, `vec` has wrong entries (either really small or really large) while with the v4.8 I obtain the correct eigenvector shape.

Best,

Lucas

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 16, 2021, 5:52pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/18 "2021-03-16T17:52:57Z")

</div>

I’m actually not quite sure how that code ever worked… I’ll fix a push after diner, I’ve got it locally but I need to figure out what was previously going on.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [March 16, 2021, 7:51pm UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/19 "2021-03-16T19:51:17Z")

</div>

This should be fixed in [this](https://github.com/FreeFem/FreeFem-sources/commit/31a36ab15fdbb34942fd5bedb4becd0db7b2f616) commit, don’t really know what I was thinking of when writing that in the first place…

---

<div class="post-metadata">

### Author: ![LucasPa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/lucaspa/32/347_2.png) [@LucasPa](https://community.freefem.org/u/LucasPa)
#### Post date: [March 17, 2021, 8:53am UTC](https://community.freefem.org/t/extract-of-left-and-right-vectors-of-multidimensional-svd/844/20 "2021-03-17T08:53:16Z")

</div>

That works all fine, thanks!
