[Scilab-users] Evaluating a polynomial on a square matrix

classic Classic list List threaded Threaded
7 messages Options
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Evaluating a polynomial on a square matrix


Dear all,

Is there some way of evaluating a polynomial on a square matrix in a matrix-wise (not component-wise) fashion?

Something like horner but matrix-wise.

Thanks!

Regards,

Federico Miyara

Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: Evaluating a polynomial on a square matrix

Hello Frederico,

> De : Federico Miyara
> Envoyé : mardi 24 septembre 2019 00:25
>
> Is there some way of evaluating a polynomial on a square matrix in a
> matrix-wise (not component-wise) fashion?

You can transform your polynomial into a usual external function Then apply it to a matrix.

e.g.

\\ **********

P = %s^2 + 2*%s + 3 \\ example of polynomial

C = coeff(P)

n = size(C, "c")

stringP = "Y = "+string(C(1))

for i = 2:n
    stringP = stringP+" + "+string(C(i))+"*x^"+string(i-1) end

deff("[Y] = fctP(x)", stringP) \\ creation of the function corresponding to the polynomial

\\ test

A = [1 2 ; 3 4]

fctP(A)

\\ [12.   17. ;  24.   33.]

\\ **********

General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: Evaluating a polynomial on a square matrix

In reply to this post by fmiyara
I'm aware that this is an inefficient way to compute the result
but unfortunately, the Horner form of polynomial does not work on matrices:

(A + 2)*A is not the same as A*A + 2*A for Scilab.

It is possible to write a better function manually,
by first calculating the successive powers before applying coefficients and adding them.

It is probably possible to create such functions automatically
by concatenating strings but it is not straightforward.

HTH,

regards

> stringP = "Y = "+string(C(1))
>
> for i = 2:n
>     stringP = stringP+" + "+string(C(i))+"*x^"+string(i-1) end
>
> deff("[Y] = fctP(x)", stringP) \\ creation of the function corresponding to the
> polynomial

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: Evaluating a polynomial on a square matrix

Hi,

I've got some time this morning (-: so this one seems to work and should be a bit more efficient:

// **********

P = %s^2 + 2*%s + 3

C = coeff(P)

n = size(C, "c")

if n > 1 then

    stringP = "X1 = x; "

    for i = 3:n

        stringP = stringP+"X"+string(i-1)+" = X"+string(i-2)+"*x; "

    end

else

    stringP = ""

end

stringP = stringP+"Y = "+string(C(1))

for i = 2:n

    stringP = stringP+" + "+string(C(i))+"*X"+string(i-1)

end

deff("[Y] = fctP(x)", stringP)

// **********

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer


General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: Evaluating a polynomial on a square matrix

Why just not writing Horner's algorithm directly ?

function out=horner_mat(p, X)
    a = coeff(p)
    out = zeros(X);
    for k = degree(p):-1:0
        out = out*X + a(k+1)*eye()
    end
endfunction

--> a=[1 2;3 4]
 a  = 

   1.   2.
   3.   4.


--> horner_mat(1+%s-%s^2,a)
 ans  =

  -5.   -8. 
  -12.  -17.


--> eye() + a - a^2
 ans  =

  -5.   -8. 
  -12.  -17.

S.
Le 24/09/2019 à 10:31, Dang Ngoc Chan, Christophe a écrit :
Hi,

I've got some time this morning (-: so this one seems to work and should be a bit more efficient:

// **********

P = %s^2 + 2*%s + 3

C = coeff(P)

n = size(C, "c")

if n > 1 then

    stringP = "X1 = x; "

    for i = 3:n

        stringP = stringP+"X"+string(i-1)+" = X"+string(i-2)+"*x; "

    end

else

    stringP = ""

end

stringP = stringP+"Y = "+string(C(1))

for i = 2:n

    stringP = stringP+" + "+string(C(i))+"*X"+string(i-1)

end

deff("[Y] = fctP(x)", stringP)

// **********

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer


General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

Re: Evaluating a polynomial on a square matrix


Stéphane,

Thanks, this is great, the use of eye is the way to go to overcome the problem of an earlier suggestion which discarded, for instance, A + 2 for not being what Horner expects; but A + 2*eye(A) is!

Thanks also for making me aware of the simplified syntax 2*eye() which adapts its size to the matrix it is added to, and is also computationally more efficient, important especially for very large matrices.

Best regards,

Federico Miyara


On 24/09/2019 05:34, Stéphane Mottelet wrote:

Why just not writing Horner's algorithm directly ?

function out=horner_mat(p, X)
    a = coeff(p)
    out = zeros(X);
    for k = degree(p):-1:0
        out = out*X + a(k+1)*eye()
    end
endfunction

--> a=[1 2;3 4]
 a  = 

   1.   2.
   3.   4.


--> horner_mat(1+%s-%s^2,a)
 ans  =

  -5.   -8. 
  -12.  -17.


--> eye() + a - a^2
 ans  =

  -5.   -8. 
  -12.  -17.

S.
Le 24/09/2019 à 10:31, Dang Ngoc Chan, Christophe a écrit :
Hi,

I've got some time this morning (-: so this one seems to work and should be a bit more efficient:

// **********

P = %s^2 + 2*%s + 3

C = coeff(P)

n = size(C, "c")

if n > 1 then

    stringP = "X1 = x; "

    for i = 3:n

        stringP = stringP+"X"+string(i-1)+" = X"+string(i-2)+"*x; "

    end

else

    stringP = ""

end

stringP = stringP+"Y = "+string(C(1))

for i = 2:n

    stringP = stringP+" + "+string(C(i))+"*X"+string(i-1)

end

deff("[Y] = fctP(x)", stringP)

// **********

Regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer


General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Christophe Dang Ngoc Chan Christophe Dang Ngoc Chan
Reply | Threaded
Open this post in threaded view
|

Re: {EXT} Re: Evaluating a polynomial on a square matrix

Hello,

> De : Federico Miyara
> Envoyé : mardi 24 septembre 2019 17:09
>
> Stéphane,
> Thanks, this is great

Indeed, my own solution proves to be poor because the conversion to a string causes a drastic truncation.

A quick example with something like P = poly(0.01:0.01:1, "s") shows a great difference with horner() (e.g. with the 1 × 1 matrix [1]).

Thanx and regards

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users