Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #60572 > unrolled thread

fast matrix determinant

Back to article view | Back to comp.lang.forth

Contents

  fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-10 10:29 -0700
    Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 05:41 -0700
      Re: fast matrix determinant mhx@iae.nl - 2017-06-12 06:03 -0700
        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 07:14 -0700
          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 07:28 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-12 11:21 -0700
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-12 16:39 -0700
          Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-14 18:05 +0300
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-14 11:05 -0700
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-14 21:37 +0300
            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 09:59 +0300
              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 03:30 -0700
                Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 17:57 +0300
                  Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 19:06 -0700
                    Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 10:21 +0300
                      Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 03:38 -0700
                        Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 15:58 +0300
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 06:04 -0700
                            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 16:19 +0300
                      Re: fast matrix determinant mhx@iae.nl - 2017-06-16 03:58 -0700
                        Re: fast matrix determinant albert@cherry.spenarnc.xs4all.nl (Albert van der Horst) - 2017-06-16 13:39 +0200
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:58 -0700
                        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:12 -0700
                          Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 05:40 -0700
                            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 17:50 +0300
                              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-16 13:39 -0700
                                Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-17 10:01 +0300
                            Re: fast matrix determinant mhx@iae.nl - 2017-06-17 01:31 -0700
                              Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-17 15:24 -0700
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-15 17:59 +0300
            Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-16 10:11 +0300
              Re: fast matrix determinant hum9 <obi65537@gmail.com> - 2017-06-17 09:20 +0300
        Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-13 16:43 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-13 22:04 -0700
          Re: fast matrix determinant mhx@iae.nl - 2017-06-14 12:23 -0700
            Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-14 18:31 -0700
      Re: fast matrix determinant krishna.myneni@ccreweb.org - 2017-06-15 19:13 -0700

Page 2 of 2 — ← Prev page 1 [2]

#60812

In article <03614458-497e-4859-b8e0-89787dfbd353@googlegroups.com>,
 <mhx@iae.nl> wrote:
>On Friday, June 16, 2017 at 9:21:21 AM UTC+2, hum9 wrote:
>> On 16.06.2017 05:06, krishna.myneni@ccreweb.org wrote:
>> > On Thursday, June 15, 2017 at 9:57:56 AM UTC-5, hum9 wrote:
>> >> On 15.06.2017 13:30, krishna...org wrote:
>...
>> I reformulated that theorem so it looks now:
>>
>> If real NxN matrix [A] has N orthonormal eigenvectors,
>> with matrix [C] of orthonormal eigenvectors as columns.
>>
>> In the case of small perturbation of [A], for every
>> unperturbed eigenvector {x} corresponds an perturbed
>> eigenvector {x+dx} collinear with {x}.
>>
>> Changing basis by [C] to a small perturbed version of [A]
>> will result an slight perturbed diagonal matrix with small
>> perturbed eigenvalues on its diagonal.
>
>I don't understand your formulation, sorry.
>A = X D X^-1, where X is matrix of eigenvectors and D the
>diagonal matrix of A's eigenvalues.
>Isn't it obvious that a (very) small perturbation of A
>will result in (very) small perturbations of D and X?
>I don't even need orthonormality for that.
>
>Interpreted physically, in my domain, it means that if
>the system matrix of a 'reasonable' application (like
>an electronic low-pass filter) changes slightly, the
>eigenvectors and eigenvalues won't change dramatically.
>I would even expect collinearity of the eigenvectors.
>
>The theorem would be more useful if you could extend it
>and say something about the relationship between small
>changes in A and the resulting changes in D and X.

It seems to me that you guys are looking for the
Gershgorin circle theorem

https://en.wikipedia.org/wiki/Gershgorin_circle_theorem

This quantifies how for off the eigenvalues can go by adding
off diagonal elements.

Not that I understood much about the discussion.

Your intuition about the effect of small perturburations is right
but hard to prove. The Gershgorin theorem suggests that
things get mixed up as soon as eigenvalues are equal or almost
equal.

>
>-marcel
-- 
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

[toc] | [prev] | [next] | [standalone]

#60819

On Friday, June 16, 2017 at 6:38:57 AM UTC-5, Albert van der Horst wrote:
...
> It seems to me that you guys are looking for the
> Gershgorin circle theorem
> 
> https://en.wikipedia.org/wiki/Gershgorin_circle_theorem
> 
> This quantifies how for off the eigenvalues can go by adding
> off diagonal elements.
> 
> Not that I understood much about the discussion.
> 
> Your intuition about the effect of small perturburations is right
> but hard to prove. ...

Ok, I realize now that what hum9 is demonstrating is
just first-order perturbation theory. The eigenvectors
are assumed to be stationary, and used to compute the
change in eigenvalues. It was difficult for me to
recognize this when presented in matrix form -- in QM
we are used to writing 1st order perturbation theory
with eigenfunctions.

Krishna

[toc] | [prev] | [next] | [standalone]

#60815

On Friday, June 16, 2017 at 5:59:01 AM UTC-5, m...@iae.nl wrote:
...
> The theorem would be more useful if you could extend it 
> and say something about the relationship between small 
> changes in A and the resulting changes in D and X.
> 

Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B

  B = C^T dA C 

provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.

Krishna

[toc] | [prev] | [next] | [standalone]

#60817

On Friday, June 16, 2017 at 7:12:34 AM UTC-5, krishna...org wrote:
> On Friday, June 16, 2017 at 5:59:01 AM UTC-5, m...@iae.nl wrote:
> ...
> > The theorem would be more useful if you could extend it 
> > and say something about the relationship between small 
> > changes in A and the resulting changes in D and X.
> > 
> 
> Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B
> 
>   B = C^T dA C 
> 
> provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.


A clearer way to see this, and which also makes it
obvious, is to write

    D  =  C^T A C

where C is the matrix of eigenvectors.

    D + dD ≈ C^T A + dA C

Thus,

    dD  ≈  C^T dA C

where only the diagonal elements are retained for the resulting matrix on the R.H.S.

In this form, the result can obviously be extended to real matrices which do not have orthogonal eigenvectors, by writing,

     dD ≈ C^-1 dA C

Krishna

[toc] | [prev] | [next] | [standalone]

#60832

On 16.06.2017 15:40, krishna.myneni@ccreweb.org wrote:
> On Friday, June 16, 2017 at 7:12:34 AM UTC-5, krishna...org wrote:
>> On Friday, June 16, 2017 at 5:59:01 AM UTC-5, m...@iae.nl wrote:
>> ...
>>> The theorem would be more useful if you could extend it
>>> and say something about the relationship between small
>>> changes in A and the resulting changes in D and X.
>>>
>>
>> Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B
>>
>>    B = C^T dA C
>>
>> provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.
> 
> 
> A clearer way to see this, and which also makes it
> obvious, is to write
> 
>      D  =  C^T A C
> 
> where C is the matrix of eigenvectors.
> 
>      D + dD ≈ C^T A + dA C
> 
> Thus,
> 
>      dD  ≈  C^T dA C
> 
> where only the diagonal elements are retained for the resulting matrix on the R.H.S.
> 
> In this form, the result can obviously be extended to real matrices which do not have orthogonal eigenvectors, by writing,
> 
>       dD ≈ C^-1 dA C
> 
> Krishna
> 
> 
> 

In mentioned context, yes, that is a *quantitative* one,
plus another *qualitative* consequence:
for every eigenvector of unperturbed matrix [A],
difference between unperturbed eigenvector {x} and
perturbed eigenvector {x+dx} is only by magnitude,
not by orientation.

hum9

[toc] | [prev] | [next] | [standalone]

#60854

On Friday, June 16, 2017 at 9:50:21 AM UTC-5, hum9 wrote:
> On 16.06.2017 15:40, krishna...org wrote:
...
> >> Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B
> >>
> >>    B = C^T dA C
> >>
> >> provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.
> > 
> > 
> > A clearer way to see this, and which also makes it
> > obvious, is to write
> > 
> >      D  =  C^T A C
> > 
> > where C is the matrix of eigenvectors.
> > 
> >      D + dD ≈ C^T A + dA C
> > 
> > Thus,
> > 
> >      dD  ≈  C^T dA C
> > 
> > where only the diagonal elements are retained for the resulting matrix on the R.H.S.
> > 
> > In this form, the result can obviously be extended to real matrices which do not have orthogonal eigenvectors, by writing,
> > 
> >       dD ≈ C^-1 dA C
> > 
> > Krishna
> > 
> > 
> > 
> 
> In mentioned context, yes, that is a *quantitative* one,
> plus another *qualitative* consequence:
> for every eigenvector of unperturbed matrix [A],
> difference between unperturbed eigenvector {x} and
> perturbed eigenvector {x+dx} is only by magnitude,
> not by orientation.
> 

Your "qualitative" result may be true in first-order, for the specific matrix case considered, but it will not be true for general perturbations for 2nd order and higher order. I guess I don't understand what you consider to be the essential result. I had interpreted the essential result to be the first-order approximation to the eigenvalues of the perturbed matrix. This is a result from perturbation theory, which goes back a long way. The wikipedia page on eigenvalue perturbation theory has some issues, but you can find your result on that page:

https://en.wikipedia.org/wiki/Eigenvalue_perturbation

The general theory provides a means of calculating not only the first order eigenvalues of the perturbed matrix, but also the first order eigenvectors of the perturbed matrix. An iterative procedure can be used to compute higher order corrections. If I find a decent reference online to point you to, I'll post it.

Regards,
Krishna

[toc] | [prev] | [next] | [standalone]

#60886

On 16.06.2017 23:39, krishna.myneni@ccreweb.org wrote:
> On Friday, June 16, 2017 at 9:50:21 AM UTC-5, hum9 wrote:
>> On 16.06.2017 15:40, krishna...org wrote:
> ...
>>>> Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B
>>>>
>>>>     B = C^T dA C
>>>>
>>>> provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.
>>>
>>>
>>> A clearer way to see this, and which also makes it
>>> obvious, is to write
>>>
>>>       D  =  C^T A C
>>>
>>> where C is the matrix of eigenvectors.
>>>
>>>       D + dD ≈ C^T A + dA C
>>>
>>> Thus,
>>>
>>>       dD  ≈  C^T dA C
>>>
>>> where only the diagonal elements are retained for the resulting matrix on the R.H.S.
>>>
>>> In this form, the result can obviously be extended to real matrices which do not have orthogonal eigenvectors, by writing,
>>>
>>>        dD ≈ C^-1 dA C
>>>
>>> Krishna
>>>
>>>
>>>
>>
>> In mentioned context, yes, that is a *quantitative* one,
>> plus another *qualitative* consequence:
>> for every eigenvector of unperturbed matrix [A],
>> difference between unperturbed eigenvector {x} and
>> perturbed eigenvector {x+dx} is only by magnitude,
>> not by orientation.
>>
> 
> Your "qualitative" result may be true in first-order, for the specific matrix case considered, but it will not be true for general perturbations for 2nd order and higher order. I guess I don't understand what you consider to be the essential result. I had interpreted the essential result to be the first-order approximation to the eigenvalues of the perturbed matrix. This is a result from perturbation theory, which goes back a long way. The wikipedia page on eigenvalue perturbation theory has some issues, but you can find your result on that page:
> 
> https://en.wikipedia.org/wiki/Eigenvalue_perturbation
> 
> The general theory provides a means of calculating not only the first order eigenvalues of the perturbed matrix, but also the first order eigenvectors of the perturbed matrix. An iterative procedure can be used to compute higher order corrections. If I find a decent reference online to point you to, I'll post it.
> 
> Regards,
> Krishna
> 
> 
> 

I reviewed proof paper and found that equation
{x}k^T * ([A] -Lk*[I]) * {dx}k = 0
   I copied wrongly from paper to computer, as
([A] -Lk*[I]) * {dx}k = 0

As a "GIGO theorem" applies, collinearity result is G(arbage).
*My apologies*.


Yet is something interesting in solving for delta eigenvector {dx}k,

([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k

as Determinant([A] -Lk*[I]) = 0.

I note [B] = [C]^T * [dA] * [C] and
transformed {dx}k into orthonormal basis -> {dy}k = [C]^T*{dx}k
with [C]^-1=[C]^T.

Transforming equation in orthonormal basis of eigenvectors, give
the solution:

dYik = |{x}k| * Bik/(Lk - Li) , for Li <> Lk

undetermined for rest, 0 * dYkk = 0.

That means components perpendicular to {x}k could be calculated,
but along the {x}k could be every finite number.

I checked Wikipedia, I found that is presented a more general case,
but in general they took the same approach as I, without factoring
terms, which I did.

Thank you for kind support,
hum9

[toc] | [prev] | [next] | [standalone]

#60892

On Friday, June 16, 2017 at 2:40:10 PM UTC+2, krishna...@ccreweb.org wrote:
> On Friday, June 16, 2017 at 7:12:34 AM UTC-5, krishna...org wrote:
> > On Friday, June 16, 2017 at 5:59:01 AM UTC-5, m...@iae.nl wrote:
> > ...
> > > The theorem would be more useful if you could extend it 
> > > and say something about the relationship between small 
> > > changes in A and the resulting changes in D and X.
> > > 
> > 
> > Hum9 omitted the main result in his reformulation statement, namely, the diagonal elements of B
> > 
> >   B = C^T dA C 
> > 
> > provide the quantitative relationship between small changes in A and the changes in D, to first order. This appears to be rigorous (and provable) when the eigenvectors of A are orthogonal.
> 
> 
> A clearer way to see this, and which also makes it
> obvious, is to write
> 
>     D  =  C^T A C
> 
> where C is the matrix of eigenvectors.
> 
>     D + dD ≈ C^T A + dA C
> 
> Thus,
> 
>     dD  ≈  C^T dA C
> 
> where only the diagonal elements are retained for the resulting matrix on the R.H.S.
> 
> In this form, the result can obviously be extended to real matrices which do not have orthogonal eigenvectors, by writing,
> 
>      dD ≈ C^-1 dA C
> 

Here is a numerical example for an electronic circuit 
(Buck converter). I'll use Matlab.

function x = theorem(dR)

% Schematic
%     S                 Vbus
% .---/---.--L--rL---.---.
% |       |           |  |
% O Vcap  - D      C === R
% |       ^           |  |
% |       |           |  |
% =       =           =  =

%% Component values
L=50e-6; rL=0.1; C=10e-6; R=15;

%% Assume this buck converter can be described with
% x = [ iCoil    vBus     vCap ]';
  A = [ -(rL)/L  -1/L      1/L
         1/C    -1/(R*C)    0
          0        0        0  ];

%% Eigenvectors:
[X,D]=eig(A);

%% Theorem:  dD ~ X^-1 dA X
% Let's increase rL from 0.1 to 0.2

  A0 = [    -rL/L      -1/L      1/L
             1/C     -1/(R*C)     0
              0          0        0  ];
           
  A1 = [ -(rL+dR)/L   -1/L      1/L
             1/C     -1/(R*C)     0  
              0          0        0  ];

dA = A1-A0;
dD = X^-1 * dA * X;

fprintf('converter eigenvalues after perturbing rL with %g\n',dR);
display(diag(D+dD,0)');
 
[~,Dn]=eig(A1); newD=diag(Dn,0)';
fprintf('actual eigenvalues\n');
display(newD);
 
%EOF

>> theorem(0)
converter eigenvalues after perturbing rL with 0
   1.0e+04 *
  -0.4333 - 4.4660i  -0.4333 + 4.4660i   0.0000 + 0.0000i
actual eigenvalues
newD =
   1.0e+04 *
  -0.4333 - 4.4660i  -0.4333 + 4.4660i   0.0000 + 0.0000i

>> theorem(0.1)
converter eigenvalues after perturbing rL with 0.1
   1.0e+04 *
  -0.5333 - 4.4713i  -0.5333 + 4.4713i   0.0000 + 0.0000i
actual eigenvalues
newD =
   1.0e+04 *
  -0.5333 - 4.4701i  -0.5333 + 4.4701i   0.0000 + 0.0000i

>> theorem(4)
converter eigenvalues after perturbing rL with 4
   1.0e+04 *
  -4.4333 - 4.6750i  -4.4333 + 4.6750i   0.0000 + 0.0000i
actual eigenvalues
newD =
   1.0e+04 *
  -4.4333 - 2.4109i  -4.4333 + 2.4109i   0.0000 + 0.0000i

>> theorem(4.706)
converter eigenvalues after perturbing rL with 4.706
   1.0e+04 *
  -5.1393 - 4.7119i  -5.1393 + 4.7119i   0.0000 + 0.0000i
actual eigenvalues
newD =
   1.0e+04 *
   -5.2082   -5.0704         0

This works surprising well (in this case), even
for huge variations in the inductor losses.

However, if one already has the eigenvalues and
eigenvectors it is difficult to see that an
approximation to the eigenvalues has practical
value. Maybe when their calculation is very time
consuming (e.g. when using SPICE).

-marcel

[toc] | [prev] | [next] | [standalone]

#60926

On Saturday, June 17, 2017 at 3:31:56 AM UTC-5, m...@iae.nl wrote:
...
> 
> Here is a numerical example for an electronic circuit 
> (Buck converter). I'll use Matlab.
> 
> function x = theorem(dR)
> 
> % Schematic
> %     S                 Vbus
> % .---/---.--L--rL---.---.
> % |       |           |  |
> % O Vcap  - D      C === R
> % |       ^           |  |
> % |       |           |  |
> % =       =           =  =
> 
> %% Component values
> L=50e-6; rL=0.1; C=10e-6; R=15;
> 
> %% Assume this buck converter can be described with
> % x = [ iCoil    vBus     vCap ]';
>   A = [ -(rL)/L  -1/L      1/L
>          1/C    -1/(R*C)    0
>           0        0        0  ];
> 
> %% Eigenvectors:
> [X,D]=eig(A);
> 
> %% Theorem:  dD ~ X^-1 dA X
> % Let's increase rL from 0.1 to 0.2
> 
...
> This works surprising well (in this case), even
> for huge variations in the inductor losses.
> 
> However, if one already has the eigenvalues and
> eigenvectors it is difficult to see that an
> approximation to the eigenvalues has practical
> value. Maybe when their calculation is very time
> consuming (e.g. when using SPICE).
...

Yes, for large matrices which  need to be swept over a small range of parameters, I expect it's much faster to just do two matrix multiplications, rather than solving the eigenvalue problem for each matrix. 

In QM, first and second order perturbation theory is used frequently to obtain analytic expressions for the eigenvalues as a function of the parameters, for example, the change in energy levels of an atom in response to an external applied electric field.

Krishna

[toc] | [prev] | [next] | [standalone]

#60781

On 15.06.2017 09:59, hum9 wrote:
> On 14.06.2017 18:05, hum9 wrote:
>> On 12.06.2017 17:14, krishna.myneni@ccreweb.org wrote:
>>> Hi Marcel,
>>>
>>> I need a fast determinant to write some pedagogical code
>>> for illustrating a molecular physics calculation. It is
>>> to be used for solving the eigenvalue problem,
>>>
>>>      D( H - lambda×S ) = 0
>>>
>>> where H and S are N×N matrices, and lambda(s) are the
>>> eigenvalue(s) to be found. D() is, of course, the
>>> determinant operator. This form of the eigenvalue
>>> problem is typical in molecular physics calculations,
>>> when the basis set of eigenfunctions are non-orthogonal.
>>>
>>> The matrix S, consists of the inner products, of the N
>>> eigenfunctions, and it would simply be the identity
>>> matrix if the N eigenfunctions were orthonormal.
>>> Typically, the hard part in setting up the problem is
>>> evaluating the elements of H and S, each of the elements
>>> requiring the evaluation of a multi-dimensional integral.
>>> Nevertheless, when N becomes large, efficient evaluation
>>> of the determinants is needed, particularly in searching
>>> for the eigenvalues.
>>>
>>> A simple way of solving the equation, once the matrix
>>> elements of H and S have been computed in the basis set,
>>> is to simply use the determinant equation within a root
>>> finder, to search for the values of lambda which satisfy
>>> the equation. This is conceptually simple and straight-
>>> forward to implement, but it's probably not the most
>>> computationally efficient method -- there are probably
>>> much better techniques available in a linear algebra
>>> library. But, for pedagogical use, it's nice. Typically,
>>> we are interested in only the lowest eigenvalue.
>>>
>>> Krishna
>>>
>>
>> Recently, after your post, I have found something that
>> could be useful:
>>
>>
>> A symbol between '[' ']' represents a matrix,
>> between '{' '}' represents vector,
>> ^T   means transpose
>> ^-1  means inverse of matrix.
>>
>> Theorem (proposition):
>>
>> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
>> and [C] its transformation matrix into orthonormal basis of [A]
>> eigenvectors.
>>
>> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
>> will result delta eigenvalues dL, for every index k, 1 <= k <= N
>>
>> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
>>
>>
>>
>> <Proof/
>>
>> [A] is a symmetric matrix NxN, with N distinct eigenvalues.
>> That means:
>>
>>       1) eigenvectors {x}i, 1 <= i <= N are linear independents
>>
>>       2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
>>
>> ---
>> 3)   [A]*{x}=L*(x).
>>
>> For a very small perturbation [dA] will follows:
>>
>> [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
>>
>> [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>      =
>> L*{x} + L*{dx} + dL*{x} + dL*{dx}
>>
>> using 3)   =>
>>
>> [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>       =
>> L*{dx} + dL*{x} + dL*{dx}
>>
>> For Norm([dA]) << Norm([A]) , second order terms are negligible, so
>>
>> [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
>>
>> ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
>>
>> For every 1 <= k <= N:
>>
>> 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
>>
>> Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
>>
>> {dx}k  =  SUM(Fik * {x}i) =>
>>
>> ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
>>       =
>> SUM(Fik * ([A] - Lk*[I])*{x}i ) =
>> SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
>> SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
>> SUM(Fik * (Li - Lk)*{x}i )
>>
>> So 4)   becomes
>>
>> SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
>>
>> multiplying at left side of equality with Transpose of {x}k results
>>
>> {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * 
>> {x}k   <=>
>>
>> SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>
>>  From 2) results that for i != k  {x}k^T * {x}i = 0,
>> also for i=k,  Li - Lk = 0, so left side of equality becames 0:
>>
>> 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>
>> if we note unity eigenvector {v}k = {x}k/Norm({x}k):
>>
>> 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
>>
>>
>> Transforming bilinear form into orthonormal basis of eigenvectors,
>> equality becames:
>>
>> 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
>>
>>     ,with {e}k versor of orthonormal eigenvalues basis.
>>
>> Let's note:
>>
>> [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
>>
>> [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
>>
>> [D] = [I]*dLk - [C]^T*[dA]*[C]
>>
>>
>> Let's note  [B] = [C]^T*[dA]*[C]  ,so
>>
>>
>> [D] = [I]*dLk - [B]  ,so 5)   becames:
>>
>>
>> 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
>>
>> 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
>>
>> 0= dLk - Bkk <=>
>>
>> *************
>> * dLk = Bkk *
>> *************
>> /Proof>
>>
>> Quickly  tested for matrix
>>
>> 0 1
>> 1 0
>>
>> with small perturbation [dA]
>>
>> a11 a12
>> a21 a22
>>
>> ,aij << 1
>>
>> and it holds :)
>>
>>
>> Have a nice day,
>> hum9
> 
> I reformulated this proposed theorem with less constraints:
> 
> If real NxN matrix [A] has N orthonormal eigenvectors,
> for a small perturbation [dA] corresponds small eigenvalues changes:
>      dLi = Bii , where [B] = [C]^T * [dA] * [C],
> with matrix [C] of orthonormal eigenvectors as columns.
> 
> hum9

Reformulated again:

If real NxN matrix [A] has N orthonormal eigenvectors,
with matrix [C] of orthonormal eigenvectors as columns.

Changing basis by [C] to a small perturbed version of [A]
will result an slight perturbed diagonal matrix with small
perturbed eigenvalues on its diagonal.

[toc] | [prev] | [next] | [standalone]

#60801

On 14.06.2017 18:05, hum9 wrote:
> On 12.06.2017 17:14, krishna.myneni@ccreweb.org wrote:
>> Hi Marcel,
>>
>> I need a fast determinant to write some pedagogical code
>> for illustrating a molecular physics calculation. It is
>> to be used for solving the eigenvalue problem,
>>
>>      D( H - lambda×S ) = 0
>>
>> where H and S are N×N matrices, and lambda(s) are the
>> eigenvalue(s) to be found. D() is, of course, the
>> determinant operator. This form of the eigenvalue
>> problem is typical in molecular physics calculations,
>> when the basis set of eigenfunctions are non-orthogonal.
>>
>> The matrix S, consists of the inner products, of the N
>> eigenfunctions, and it would simply be the identity
>> matrix if the N eigenfunctions were orthonormal.
>> Typically, the hard part in setting up the problem is
>> evaluating the elements of H and S, each of the elements
>> requiring the evaluation of a multi-dimensional integral.
>> Nevertheless, when N becomes large, efficient evaluation
>> of the determinants is needed, particularly in searching
>> for the eigenvalues.
>>
>> A simple way of solving the equation, once the matrix
>> elements of H and S have been computed in the basis set,
>> is to simply use the determinant equation within a root
>> finder, to search for the values of lambda which satisfy
>> the equation. This is conceptually simple and straight-
>> forward to implement, but it's probably not the most
>> computationally efficient method -- there are probably
>> much better techniques available in a linear algebra
>> library. But, for pedagogical use, it's nice. Typically,
>> we are interested in only the lowest eigenvalue.
>>
>> Krishna
>>
> 
> Recently, after your post, I have found something that
> could be useful:
> 
> 
> A symbol between '[' ']' represents a matrix,
> between '{' '}' represents vector,
> ^T   means transpose
> ^-1  means inverse of matrix.
> 
> Theorem (proposition):
> 
> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
> and [C] its transformation matrix into orthonormal basis of [A]
> eigenvectors.
> 
> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
> will result delta eigenvalues dL, for every index k, 1 <= k <= N
> 
> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
> 
> 
> 
> <Proof/
> 
> [A] is a symmetric matrix NxN, with N distinct eigenvalues.
> That means:
> 
>       1) eigenvectors {x}i, 1 <= i <= N are linear independents
> 
>       2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
> 
> ---
> 3)   [A]*{x}=L*(x).
> 
> For a very small perturbation [dA] will follows:
> 
> [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
> 
> [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
>      =
> L*{x} + L*{dx} + dL*{x} + dL*{dx}
> 
> using 3)   =>
> 
> [A]*{dx} + [dA]*{x} + [dA]*{dx}
>       =
> L*{dx} + dL*{x} + dL*{dx}
> 
> For Norm([dA]) << Norm([A]) , second order terms are negligible, so
> 
> [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
> 
> ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
> 
> For every 1 <= k <= N:
> 
> 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
> 
> Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
> 
> {dx}k  =  SUM(Fik * {x}i) =>
> 
> ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
>       =
> SUM(Fik * ([A] - Lk*[I])*{x}i ) =
> SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
> SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
> SUM(Fik * (Li - Lk)*{x}i )
> 
> So 4)   becomes
> 
> SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
> 
> multiplying at left side of equality with Transpose of {x}k results
> 
> {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * {x}k   <=>
> 
> SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
> 
>  From 2) results that for i != k  {x}k^T * {x}i = 0,
> also for i=k,  Li - Lk = 0, so left side of equality becames 0:
> 
> 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
> 
> if we note unity eigenvector {v}k = {x}k/Norm({x}k):
> 
> 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
> 
> 
> Transforming bilinear form into orthonormal basis of eigenvectors,
> equality becames:
> 
> 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
> 
>     ,with {e}k versor of orthonormal eigenvalues basis.
> 
> Let's note:
> 
> [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
> 
> [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
> 
> [D] = [I]*dLk - [C]^T*[dA]*[C]
> 
> 
> Let's note  [B] = [C]^T*[dA]*[C]  ,so
> 
> 
> [D] = [I]*dLk - [B]  ,so 5)   becames:
> 
> 
> 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
> 
> 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
> 
> 0= dLk - Bkk <=>
> 
> *************
> * dLk = Bkk *
> *************
> /Proof>
> 
> Quickly  tested for matrix
> 
> 0 1
> 1 0
> 
> with small perturbation [dA]
> 
> a11 a12
> a21 a22
> 
> ,aij << 1
> 
> and it holds :)
> 
> 
> Have a nice day,
> hum9


It was demonstrated that left side of equality 4) is zero:
	([A] -Lk*[I]) * {dx}k = 0
That means delta-eigenvector is also an eigenvector of [A]!

That means for an eigenvalue Lk of [A], perturbed eigenvector
{x+dx}k is collinear with unperturbed eigenvector of {x}k.

So theorem can be reformulated:

<Theorem>
If real NxN matrix [A] has N orthonormal eigenvectors,
with matrix [C] of orthonormal eigenvectors as columns.

In the case of small perturbation of [A], for every
unperturbed eigenvector {x} corresponds an perturbed
eigenvector {x+dx} collinear with {x}.

Changing basis by [C] to a small perturbed version of [A]
will result an slight perturbed diagonal matrix with small
perturbed eigenvalues on its diagonal.
</Theorem>

[toc] | [prev] | [next] | [standalone]

#60885

On 16.06.2017 10:11, hum9 wrote:
> On 14.06.2017 18:05, hum9 wrote:
>> On 12.06.2017 17:14, krishna.myneni@ccreweb.org wrote:
>>> Hi Marcel,
>>>
>>> I need a fast determinant to write some pedagogical code
>>> for illustrating a molecular physics calculation. It is
>>> to be used for solving the eigenvalue problem,
>>>
>>>      D( H - lambda×S ) = 0
>>>
>>> where H and S are N×N matrices, and lambda(s) are the
>>> eigenvalue(s) to be found. D() is, of course, the
>>> determinant operator. This form of the eigenvalue
>>> problem is typical in molecular physics calculations,
>>> when the basis set of eigenfunctions are non-orthogonal.
>>>
>>> The matrix S, consists of the inner products, of the N
>>> eigenfunctions, and it would simply be the identity
>>> matrix if the N eigenfunctions were orthonormal.
>>> Typically, the hard part in setting up the problem is
>>> evaluating the elements of H and S, each of the elements
>>> requiring the evaluation of a multi-dimensional integral.
>>> Nevertheless, when N becomes large, efficient evaluation
>>> of the determinants is needed, particularly in searching
>>> for the eigenvalues.
>>>
>>> A simple way of solving the equation, once the matrix
>>> elements of H and S have been computed in the basis set,
>>> is to simply use the determinant equation within a root
>>> finder, to search for the values of lambda which satisfy
>>> the equation. This is conceptually simple and straight-
>>> forward to implement, but it's probably not the most
>>> computationally efficient method -- there are probably
>>> much better techniques available in a linear algebra
>>> library. But, for pedagogical use, it's nice. Typically,
>>> we are interested in only the lowest eigenvalue.
>>>
>>> Krishna
>>>
>>
>> Recently, after your post, I have found something that
>> could be useful:
>>
>>
>> A symbol between '[' ']' represents a matrix,
>> between '{' '}' represents vector,
>> ^T   means transpose
>> ^-1  means inverse of matrix.
>>
>> Theorem (proposition):
>>
>> Be a real symmetric matrix [A] NxN  with N distinct nonzero eigenvalues
>> and [C] its transformation matrix into orthonormal basis of [A]
>> eigenvectors.
>>
>> Then, for very small perturbation [dA] (Norm([dA] << Norm([A]))
>> will result delta eigenvalues dL, for every index k, 1 <= k <= N
>>
>> dLk = Bkk  ,where [B] = [C]^T * [dA] * [C]
>>
>>
>>
>> <Proof/
>>
>> [A] is a symmetric matrix NxN, with N distinct eigenvalues.
>> That means:
>>
>>       1) eigenvectors {x}i, 1 <= i <= N are linear independents
>>
>>       2) for every i j,   1 <= i != j <=N,   {x}i^T * {x}j = 0.
>>
>> ---
>> 3)   [A]*{x}=L*(x).
>>
>> For a very small perturbation [dA] will follows:
>>
>> [A + dA] * {x + dx} = (L + dl) * {x + dx} <=>
>>
>> [A]*{x} + [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>      =
>> L*{x} + L*{dx} + dL*{x} + dL*{dx}
>>
>> using 3)   =>
>>
>> [A]*{dx} + [dA]*{x} + [dA]*{dx}
>>       =
>> L*{dx} + dL*{x} + dL*{dx}
>>
>> For Norm([dA]) << Norm([A]) , second order terms are negligible, so
>>
>> [A]*{dx} + [dA]*{x} = L*{dx} + dL*{x} <=>
>>
>> ([A] - L*[I])*{dx} = ([I]*dL - [dA])*{x} <=>
>>
>> For every 1 <= k <= N:
>>
>> 4)   ([A] -Lk*[I]) * {dx}k = ([I]*dLk - [dA]) * {x}k
>>
>> Using 1), vector {dx}k can be expressed in terms of eigenvectors {x}i, so
>>
>> {dx}k  =  SUM(Fik * {x}i) =>
>>
>> ([A] - Lk*[I]) * {dx}k = ([A] - Lk*[I]) * SUM(Fik * {x}i)
>>       =
>> SUM(Fik * ([A] - Lk*[I])*{x}i ) =
>> SUM(Fik * ([A]*{x}i - Lk*[I]*{x}i) ) =
>> SUM(Fik * (Li*{x}i  - Lk*{x}i) ) =
>> SUM(Fik * (Li - Lk)*{x}i )
>>
>> So 4)   becomes
>>
>> SUM(Fik * (Li - Lk)*{x}i ) = ([I]*dLk - [dA]) * {x}k
>>
>> multiplying at left side of equality with Transpose of {x}k results
>>
>> {x}k^T * SUM(Fik * (Li - Lk)*{x}i  = {x}k^T * ([I]*dLk - [dA]) * 
>> {x}k   <=>
>>
>> SUM(Fik * (Li -Lk) * {x}k^T * {x}i = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>
>>  From 2) results that for i != k  {x}k^T * {x}i = 0,
>> also for i=k,  Li - Lk = 0, so left side of equality becames 0:
>>
>> 0 = {x}k^T * ([I]*dLk - [dA]) * {x}k
>>
>> if we note unity eigenvector {v}k = {x}k/Norm({x}k):
>>
>> 0 = {v}k^T *  ([I]*dLk - [dA]) * {v}k
>>
>>
>> Transforming bilinear form into orthonormal basis of eigenvectors,
>> equality becames:
>>
>> 5)   0 = {e}k^T * ([C]^T * ([I]*dLk - [dA]) * [C]) * {e}k
>>
>>     ,with {e}k versor of orthonormal eigenvalues basis.
>>
>> Let's note:
>>
>> [D] = [C]^T * ([I]*dLk - [dA]) * [C]  <=>
>>
>> [D] = [C]^T*[I]*[C] * dLk  -  [C]^T*[dA]*[C]  <=>
>>
>> [D] = [I]*dLk - [C]^T*[dA]*[C]
>>
>>
>> Let's note  [B] = [C]^T*[dA]*[C]  ,so
>>
>>
>> [D] = [I]*dLk - [B]  ,so 5)   becames:
>>
>>
>> 0 = {e}k^T * ([I]*dLk - [B]) * {e}k  <=>
>>
>> 0 = {e}k^T * [I] * {e}k * dLk  -  {e}k^T * [B] * {e}k  <=>
>>
>> 0= dLk - Bkk <=>
>>
>> *************
>> * dLk = Bkk *
>> *************
>> /Proof>
>>
>> Quickly  tested for matrix
>>
>> 0 1
>> 1 0
>>
>> with small perturbation [dA]
>>
>> a11 a12
>> a21 a22
>>
>> ,aij << 1
>>
>> and it holds :)
>>
>>
>> Have a nice day,
>> hum9
> 
> 
> It was demonstrated that left side of equality 4) is zero:
>      ([A] -Lk*[I]) * {dx}k = 0

I reviewed proof paper and I found that:
{x}k^T * ([A] -Lk*[I]) * {dx}k = 0
which means that collinearity is not part is false.

> That means delta-eigenvector is also an eigenvector of [A]! >
> That means for an eigenvalue Lk of [A], perturbed eigenvector
> {x+dx}k is collinear with unperturbed eigenvector of {x}k.
> 
That part is false.

> So theorem can be reformulated:
> 
> <Theorem>
> If real NxN matrix [A] has N orthonormal eigenvectors,
> with matrix [C] of orthonormal eigenvectors as columns.
>

<False>


> In the case of small perturbation of [A], for every
> unperturbed eigenvector {x} corresponds an perturbed
> eigenvector {x+dx} collinear with {x}.
> 
</False>

> Changing basis by [C] to a small perturbed version of [A]
> will result an slight perturbed diagonal matrix with small
> perturbed eigenvalues on its diagonal.
> </Theorem>

[toc] | [prev] | [next] | [standalone]

#60708

On Monday, June 12, 2017 at 8:03:44 AM UTC-5, m...@iae.nl wrote:
...
> ... In fact, I have been looking for a PLACE
> (Matlab) replacement for years. SLICOT has it in
> (Fortran) source, but I didn't get to converting 
> (or linking) it yet.

So, the naive answer to implementing PLACE seems obvious, but I would guess there is a serious problem with numerical stability, particularly with respect to matrix inversion. If I understand the problem correctly, you are given A and B matrices, and you need to find a matrix K, such that A - BK has a specified set of eigenvalues {p0, p1, p2, ... }.

Let P be a diagonal matrix with elements, p0, p1, p2, ... along the diagonal, and zero everywhere else. Then,

    P = A - BK

    A - P = BK

    K = B^-1( A - P )

where B^-1 is the inverse matrix. As I recall, Numerical Recipes has a routine for dealing with nearly singular matrices, so that might be useful.

Krishna

[toc] | [prev] | [next] | [standalone]

#60710

On Wednesday, June 14, 2017 at 1:43:04 AM UTC+2, krishna...@ccreweb.org wrote:
> On Monday, June 12, 2017 at 8:03:44 AM UTC-5, m...@iae.nl wrote:
> ...
> > ... In fact, I have been looking for a PLACE
> > (Matlab) replacement for years. SLICOT has it in
> > (Fortran) source, but I didn't get to converting 
> > (or linking) it yet.
> 
> So, the naive answer to implementing PLACE seems obvious, but I would
> guess there is a serious problem with numerical stability,
> particularly with respect to matrix inversion. If I understand the
> problem correctly, you are given A and B matrices, and you need
> to find a matrix K, such that A - BK has a specified set of
> eigenvalues {p0, p1, p2, ... }.
> 
> Let P be a diagonal matrix with elements, p0, p1, p2, ... along the diagonal,
> and zero everywhere else. Then,
>
>     P = A - BK
>
>     A - P = BK
>
>     K = B^-1( A - P )
>
> where B^-1 is the inverse matrix. As I recall, Numerical Recipes has a
> routine for dealing with nearly singular matrices, so that might be useful.

The problem I had was when B has more than a single column.
In that case there are many solutions, and one needs to search for the
most robust one. I encountered problems with orthogonal projections
in the complex plane that I couldn't solve at that time.

However, I did find the Forth solution when m (columns) is 1.

DOC
(*
	evplace ( A{ B{ P{ K{ -- ) computes a state-feedback matrix K such that
	the eigenvalues of A-B*K are those specified in the vector P.

	Ref1:: Kautsky, Nichols, Van Dooren, "Robust Pole Assignment in Linear 
	       State Feedback," Intl. J. Control, 41(1985)5, pp 1129-1155
	Ref2:: Richard J. Vaccaro, "Digital Control - A State-Space Approach",
	       pp 389-398, 412-415

	The refs present code in rather dense Matlab form. No translation to
	other languages is known. 
*)
ENDDOC

-marcel

[toc] | [prev] | [next] | [standalone]

#60762

On Wednesday, June 14, 2017 at 1:43:04 AM UTC+2, krishna...@ccreweb.org wrote:
> On Monday, June 12, 2017 at 8:03:44 AM UTC-5, m...@iae.nl wrote:
> ...
> > ... In fact, I have been looking for a PLACE
> > (Matlab) replacement for years. SLICOT has it in
> > (Fortran) source, but I didn't get to converting 
> > (or linking) it yet.
> 
> So, the naive answer to implementing PLACE seems obvious,
[..]
> 
> Let P be a diagonal matrix with elements, p0, p1, p2, ... 
> along the diagonal, and zero everywhere else. Then,
> 
>     P = A - BK
> 
>     A - P = BK
> 
>     K = B^-1( A - P )
> 
> where B^-1 is the inverse matrix. 

Don't you mean  XPX^-1 = A - BK ? 
X is the matrix of eigenvectors of A - BK (that aren't given).

B is (in general) not a square matrix and thus has no 
inverse.

-marcel

[toc] | [prev] | [next] | [standalone]

#60765

On Wednesday, June 14, 2017 at 2:23:56 PM UTC-5, m...@iae.nl wrote:
> On Wednesday, June 14, 2017 at 1:43:04 AM UTC+2, krishna...org wrote:
> > On Monday, June 12, 2017 at 8:03:44 AM UTC-5, m...@iae.nl wrote:
> > ...
> > > ... In fact, I have been looking for a PLACE
> > > (Matlab) replacement for years. SLICOT has it in
> > > (Fortran) source, but I didn't get to converting 
> > > (or linking) it yet.
> > 
> > So, the naive answer to implementing PLACE seems obvious,
> [..]
> > 
> > Let P be a diagonal matrix with elements, p0, p1, p2, ... 
> > along the diagonal, and zero everywhere else. Then,
> > 
> >     P = A - BK
> > 
> >     A - P = BK
> > 
> >     K = B^-1( A - P )
> > 
> > where B^-1 is the inverse matrix. 
> 
> Don't you mean  XPX^-1 = A - BK ? 
> X is the matrix of eigenvectors of A - BK (that aren't given).
> ...

No, I meant what I said above, according to the following
logic. First, let's assume B is a square matrix of the
same size as A.

1) The matrix K is initially undetermined, and we want
   to determine K to meet the requirement that the
   matrix A - BK has the specified set of eigenvalues,
   {p0, p1, p2, ... }.

2) We are free to choose any matrix which meets the
   above requirements, so the simplest is the diagonal
   matrix P , with the specified eigenvalues along the
   diagonal, and zeros everywhere else. The matrix, P,
   is the simplest matrix which has as its eigenvalues,
   the set {p0, p1, ... }.

3) Then, setting P = A - BK, we can solve for K, as
   above. In general, this procedure will not work as
   B might be a singular matrix and not invertible.
   However, one may be able to use SVD (singular value
   decomposition) to find the pseudo-inverse of B, and
   this may provide an estimate of K.

> B is (in general) not a square matrix and thus has no 
> inverse.

Suppose A and P have sizes n x n, and B has size n x m.
Then, K has to be an m x n matrix. If we define

    K = B^T K'

where K' is an n x n matrix, then,

    K' = (B B^T)^-1 (A - P)

Again, one might have to use SVD to obtain the inverse.

I'm not at all sure that the above procedure will give a
reasonable estimate for K in your problem. But, logically,
it seems like it should work for non-singular cases.

Krishna

[toc] | [prev] | [next] | [standalone]

#60800

On Monday, June 12, 2017 at 7:41:29 AM UTC-5, krishna...org wrote:
> On Saturday, June 10, 2017 at 12:29:13 PM UTC-5, krishna...org wrote:
...
> 
> A properly factored version of kForth's fsl/extras/determ.4th will allow for optimization, both in the Forth source and by selective replacement with CODE word(s). A revised version, called determ2.4th, is given below. It is considerably faster, even though it doesn't use any CODE words, yet.
> 

The new, factored version of determ.fs, replaces the older version and may be found at

ftp://ccreweb.org/software/fsl/extras/determ.fs

The above is the version for ANS Forth with a separate fp stack. The updated kForth version is available at

ftp://ccreweb.org/software/kforth/examples/fsl/extras/determ.4th

No assembler optimizations have been performed.

Krishna

[toc] | [prev] | [standalone]

Page 2 of 2 — ← Prev page 1 [2]

Back to top | Article view | comp.lang.forth

csiph-web