[Trilinos-Users] problems with transpose, inverse and products with Epetra_CRSMatrix
Erik Boman
egboman at sandia.gov
Fri Sep 12 10:00:22 MDT 2014
Oops, I missed the important detail that you only use the diagonal of A^T A.
I still doubt if GMRES is the best way to solve least-squares problems.
You could try CG on the normal equations (CGNE or CGLS) but LSQR should
be more stable.
Erik
Heroux, Michael A wrote:
> Sorry. Yes, of course. But this is still possible. You just need to
> compute the reciprocal values on BT:
>
> BT.Reciprocal(BT);
>
> Mike
>
> On 9/12/14, 10:17 AM, "Deepak Garg" <deepak.garg at pi.ingv.it> wrote:
>
>
>> Dear Mike
>>
>> Thank you very much for your kind reply.
>>
>> I have a little doubt. We need inverse of diag(A^T A).
>>
>> B^T should be
>>
>> B^T = A * (diag(A^T A))^{-1}
>>
>> How will we compute this inverse?
>>
>> Thank you
>> Deepak
>>
>> ----- Original Message -----
>> From: "Mike Heroux" <MHeroux at CSBSJU.EDU>
>> To: "Mike Heroux" <MHeroux at CSBSJU.EDU>, "Antonella Longo"
>> <longo at pi.ingv.it>
>> Cc: trilinos-users at software.sandia.gov, "Deepak Garg"
>> <deepak.garg at pi.ingv.it>
>> Sent: Friday, September 12, 2014 4:54:18 PM
>> Subject: Re: [Trilinos-Users] problems with transpose, inverse and
>> products with Epetra_CRSMatrix
>>
>> Small error below
>>
>> sum +=curRowVals[i];
>>
>> Should be
>> sum += curRowVals[j]*curRowVals[j];
>>
>> Please note that I did not compile or test this code, so be suspicious of
>> details :).
>>
>> Mike
>>
>> On 9/12/14, 9:49 AM, "Heroux, Mike" <MHeroux at CSBSJU.EDU> wrote:
>>
>>
>>> Antonella,
>>>
>>> First, although I don’t think you will need it for this situation, you
>>> can
>>> enable EpetraExt using Cmake ENABLE option and then re-doing the make
>>> install.
>>>
>>> Second, forming B can be done much more efficiently than computing the
>>> entire normal matrix. You really just need the diagonal of the matrix,
>>> which is obtained entry by entry from the sum of squares of the matrix
>>> rows. The following code should be a good start (although I have not
>>> tested it) for computing diag(A^T A):
>>>
>>> Assuming A is your Epetra_CrsMatrix object and diagATA is an
>>> Epetra_Vector
>>> defined using the RowMap of A (Epetra_Vector diagATA(A.RowMap()))
>>>
>>> this function should compute diag(A^T A):
>>>
>>> void computeB(const Epetra_CrsMatrix & A, Epetra_Vector & diagATA) {
>>>
>>> // Preconditions: RowMap of A == DomainMap of A == RangeMap of A (These
>>> conditions are typically true)
>>>
>>> for(int i=0; i<A.NumMyRows(); i++) {
>>> int numEntries = A.Graph().NumMyIndices(i);
>>> double * curRowVals = A.Values(i);
>>> double sum = 0.0;
>>> for(int j=0; j< numEntries; j++)
>>> sum += curRowVals[j];
>>> diagATA[i] = sum;
>>> }
>>> return;
>>> }
>>>
>>> On exit, diagATA will be a vector containing the diagonal elements of A^T
>>> A.
>>>
>>>
>>> Next you can compute B^T = A * diag(A^T A) as follows:
>>>
>>> Epetra_CrsMatrix BT(A); // makes copy of A into B^T
>>> BT.RightScale(diagATA); // Scale B in-place with diag(A^T A)
>>>
>>> At this point you have the elements you need for working with the system
>>> BAx=Bb.
>>>
>>> You can continue by explicitly forming BA and Bb, if you like (not
>>> necessarily recommended, but might be OK).
>>>
>>> BA is formed by calling the MatrixMatrix functions Alicia recommended
>>> from
>>> EpetraExt.
>>>
>>> Bb is formed by:
>>>
>>> Epetra_Vector Bb(BT.DomainMap());
>>> BT.Multiply(true, b, Bb);
>>>
>>> Hopefully this is enough information to get you started.
>>>
>>> Mike
>>>
>>> On 9/11/14, 7:10 AM, "Antonella Longo" <longo at pi.ingv.it> wrote:
>>>
>>>
>>>> Thank you for your prompt answer.
>>>> In this email we try to explain the complete problem that we must solve.
>>>>
>>>> We have a non-square system of equation for which we found a
>>>> preconditioner
>>>> to make it square so that it can be solved with GMRES.
>>>>
>>>> Our initial non-square system is Ax=b.
>>>> The preconditioner requires the computation of the matrix B
>>>> that is $B=(diag(A^T A))^{-1}A^T$.
>>>> A is an Epetra_CRSMatrix and we have its map.
>>>> The final system to be solved will be BAx=Bb.
>>>>
>>>> We were not successful to make this B matrix defined above.
>>>>
>>>> Another problem that we have is that even though the EpetraExt package
>>>> is installed, and the proper ".h" files are included, the compiler does
>>>> not
>>>> recognize the functions that you suggested us.
>>>>
>>>> We are wondering if you could suggest us the right functions
>>>> to solve the above problem.
>>>>
>>>> -----------------------------------------------
>>>> Antonella Longo
>>>> Istituto Nazionale di Geofisica e Vulcanologia
>>>> Sezione di Pisa
>>>> Via Uguccione della Faggiola, 32
>>>> Pisa, Italy
>>>> Tel. office: 050 8311939
>>>> Cell.: 339 7532089
>>>>
>>>> ----- Original Message -----
>>>> From: "Mike Heroux" <MHeroux at CSBSJU.EDU>
>>>> To: "Antonella Longo" <longo at pi.ingv.it>, "Deepak Garg"
>>>> <deepak.garg at pi.ingv.it>
>>>> Cc: trilinos-users at software.sandia.gov, "Alicia Marie Klinvex"
>>>> <aklinvex at purdue.edu>
>>>> Sent: Wednesday, September 10, 2014 11:17:47 PM
>>>> Subject: Re: [Trilinos-Users] problems with transpose, inverse and
>>>> products with Epetra_CRSMatrix
>>>>
>>>> In addition to the pointers that Alicia provided, you can compute the
>>>> matrix-transpose times vector product using the specific
>>>> Epetra_CrsMatrix
>>>> method called Multiply(). This method has a bool input argument that
>>>> indicates whether or not to apply the matrix-transpose operation:
>>>>
>>>> http://trilinos.org/docs/dev/packages/epetra/doc/html/classEpetra__CrsMa
>>>> t
>>>> r
>>>> i
>>>> x.html#ae97a6c65b11aa29b3bb25d9788c56b3a
>>>>
>>>> In general, it is unusual (but not necessarily wrong) to explicitly form
>>>> the transpose of a matrix. In many situations, just being able to apply
>>>> the action of the matrix transpose is sufficient.
>>>>
>>>> Following up on Alicia¹s comments on inverses, it is almost never best
>>>> to
>>>> explicitly compute the inverse of a sparse matrix, unless it is small or
>>>> has very special structure.
>>>>
>>>> Perhaps if you describe the algorithms you are trying to implement we
>>>> can
>>>> suggest some other options.
>>>>
>>>> Mike
>>>>
>>>> On 9/10/14, 8:10 AM, "Alicia Marie Klinvex" <aklinvex at purdue.edu> wrote:
>>>>
>>>>
>>>>> Hello Antonella and Deepak,
>>>>>
>>>>> The matrix multiplication is in the EpetraExt package:
>>>>> http://trilinos.org/docs/dev/packages/epetraext/doc/html/classEpetraExt
>>>>> _
>>>>> 1
>>>>> _
>>>>> 1MatrixMatrix.html
>>>>>
>>>>> I believe the transpose function you seek is also in EpetraExt:
>>>>> http://trilinos.org/docs/dev/packages/epetraext/doc/html/classEpetraExt
>>>>> _
>>>>> 1
>>>>> _
>>>>> 1RowMatrix__Transpose.html
>>>>> (I could be wrong about that one; I've never used that particular
>>>>> function before.)
>>>>>
>>>>> The matrix vector multiplication is in Epetra (called Apply):
>>>>> http://trilinos.org/docs/dev/packages/epetra/doc/html/classEpetra__CrsM
>>>>> a
>>>>> t
>>>>> r
>>>>> ix.html#a77c4232f1b3d32ff76504ed314df5f76
>>>>>
>>>>> Do you absolutely need to compute the inverse of your matrix, or do you
>>>>> just want to apply the inverse to a vector? I'm afraid I don't know
>>>>> whether Trilinos even has a function to invert a sparse matrix, but if
>>>>> you were only interested in solving linear systems, it certainly has
>>>>> functions for that.
>>>>>
>>>>> Best wishes,
>>>>> Alicia Klinvex
>>>>>
>>>>>
>>>>>
>>>>> ----- Original Message -----
>>>>> From: "Antonella Longo" <longo at pi.ingv.it>
>>>>> To: trilinos-users at software.sandia.gov
>>>>> Cc: "Deepak Garg" <deepak.garg at pi.ingv.it>
>>>>> Sent: Wednesday, September 10, 2014 9:10:46 AM
>>>>> Subject: [Trilinos-Users] problems with transpose, inverse and products
>>>>> with Epetra_CRSMatrix
>>>>>
>>>>> We are making computations with Epetra_CRSMatrices,
>>>>> but we could not find any mathematical functions for
>>>>> these matrices.
>>>>> In particular we need to compute the transpose and
>>>>> the inverse of an Epetra_CRSMatrix, to calculate
>>>>> the product of two Epetra_CRSMatrices, and the
>>>>> product of an Epetra_CRSMatrix with an Epetra_Vector.
>>>>> Does there exist any standard functions for these operations?
>>>>> We checked the functions of the above mentioned classes but we found
>>>>> nothing.
>>>>>
>>>>> Thanks for your kind attention.
>>>>>
>>>>> Antonella Longo and Deepak Garg
>>>>>
>>>>>
>>>>>
>>>>> -----------------------------------------------
>>>>> Antonella Longo
>>>>> Istituto Nazionale di Geofisica e Vulcanologia
>>>>> Sezione di Pisa
>>>>> Via Uguccione della Faggiola, 32
>>>>> Pisa, Italy
>>>>> Tel. office: 050 8311939
>>>>> Cell.: 339 7532089
>>>>>
>>>>> _______________________________________________
>>>>> Trilinos-Users mailing list
>>>>> Trilinos-Users at software.sandia.gov
>>>>> https://software.sandia.gov/mailman/listinfo/trilinos-users
>>>>> _______________________________________________
>>>>> Trilinos-Users mailing list
>>>>> Trilinos-Users at software.sandia.gov
>>>>> https://software.sandia.gov/mailman/listinfo/trilinos-users
>>>>>
>>> _______________________________________________
>>> Trilinos-Users mailing list
>>> Trilinos-Users at software.sandia.gov
>>> https://software.sandia.gov/mailman/listinfo/trilinos-users
>>>
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> https://software.sandia.gov/mailman/listinfo/trilinos-users
>
More information about the Trilinos-Users
mailing list